Actual source code: ex7.c
petsc-3.12.2 2019-11-22
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 FormJacobianNoMatrix(SNES,Vec,Mat,Mat,void*);
9: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
10: extern PetscErrorCode FormFunctioni(void *,PetscInt,Vec,PetscScalar *);
11: extern PetscErrorCode OtherFunctionForDifferencing(void*,Vec,Vec);
12: extern PetscErrorCode FormInitialGuess(SNES,Vec);
13: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
15: typedef struct {
16: PetscViewer viewer;
17: } MonitorCtx;
19: typedef struct {
20: PetscBool variant;
21: } AppCtx;
23: int main(int argc,char **argv)
24: {
25: SNES snes; /* SNES context */
26: SNESType type = SNESNEWTONLS; /* default nonlinear solution method */
27: Vec x,r,F,U; /* vectors */
28: Mat J,B; /* Jacobian matrix-free, explicit preconditioner */
29: AppCtx user; /* user-defined work context */
30: PetscScalar h,xp = 0.0,v;
31: PetscInt its,n = 5,i;
33: PetscBool puremf = PETSC_FALSE;
35: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
36: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
37: PetscOptionsHasName(NULL,NULL,"-variant",&user.variant);
38: h = 1.0/(n-1);
40: /* Set up data structures */
41: VecCreateSeq(PETSC_COMM_SELF,n,&x);
42: PetscObjectSetName((PetscObject)x,"Approximate Solution");
43: VecDuplicate(x,&r);
44: VecDuplicate(x,&F);
45: VecDuplicate(x,&U);
46: PetscObjectSetName((PetscObject)U,"Exact Solution");
48: /* create explict matrix preconditioner */
49: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,&B);
51: /* Store right-hand-side of PDE and exact solution */
52: for (i=0; i<n; i++) {
53: v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
54: VecSetValues(F,1,&i,&v,INSERT_VALUES);
55: v = xp*xp*xp;
56: VecSetValues(U,1,&i,&v,INSERT_VALUES);
57: xp += h;
58: }
60: /* Create nonlinear solver */
61: SNESCreate(PETSC_COMM_WORLD,&snes);
62: SNESSetType(snes,type);
64: /* Set various routines and options */
65: SNESSetFunction(snes,r,FormFunction,F);
66: if (user.variant) {
67: /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */
68: MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J);
69: MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes);
70: MatMFFDSetFunctioni(J,FormFunctioni);
71: /* Use the matrix free operator for both the Jacobian used to define the linear system and used to define the preconditioner */
72: /* This tests MatGetDiagonal() for MATMFFD */
73: PetscOptionsHasName(NULL,NULL,"-puremf",&puremf);
74: } else {
75: /* create matrix free matrix for Jacobian */
76: MatCreateSNESMF(snes,&J);
77: /* demonstrates differencing a different function than FormFunction() to apply a matrix operator */
78: /* note we use the same context for this function as FormFunction, the F vector */
79: MatMFFDSetFunction(J,OtherFunctionForDifferencing,F);
80: }
82: /* Set various routines and options */
83: SNESSetJacobian(snes,J,puremf ? J : B,puremf ? FormJacobianNoMatrix : FormJacobian,&user);
84: SNESSetFromOptions(snes);
86: /* Solve nonlinear system */
87: FormInitialGuess(snes,x);
88: SNESSolve(snes,NULL,x);
89: SNESGetIterationNumber(snes,&its);
90: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
92: /* Free data structures */
93: VecDestroy(&x); VecDestroy(&r);
94: VecDestroy(&U); VecDestroy(&F);
95: MatDestroy(&J); MatDestroy(&B);
96: SNESDestroy(&snes);
97: PetscFinalize();
98: return ierr;
99: }
100: /* -------------------- Evaluate Function F(x) --------------------- */
102: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
103: {
104: const PetscScalar *xx,*FF;
105: PetscScalar *ff,d;
106: PetscInt i,n;
107: PetscErrorCode ierr;
109: VecGetArrayRead(x,&xx);
110: VecGetArray(f,&ff);
111: VecGetArrayRead((Vec) dummy,&FF);
112: VecGetSize(x,&n);
113: d = (PetscReal)(n - 1); d = d*d;
114: ff[0] = xx[0];
115: 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];
116: ff[n-1] = xx[n-1] - 1.0;
117: VecRestoreArrayRead(x,&xx);
118: VecRestoreArray(f,&ff);
119: VecRestoreArrayRead((Vec)dummy,&FF);
120: return 0;
121: }
123: PetscErrorCode FormFunctioni(void *dummy,PetscInt i,Vec x,PetscScalar *s)
124: {
125: const PetscScalar *xx,*FF;
126: PetscScalar d;
127: PetscInt n;
128: PetscErrorCode ierr;
129: SNES snes = (SNES) dummy;
130: Vec F;
132: SNESGetFunction(snes,NULL,NULL,(void**)&F);
133: VecGetArrayRead(x,&xx);
134: VecGetArrayRead(F,&FF);
135: VecGetSize(x,&n);
136: d = (PetscReal)(n - 1); d = d*d;
137: if (i == 0) {
138: *s = xx[0];
139: } else if (i == n-1) {
140: *s = xx[n-1] - 1.0;
141: } else {
142: *s = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
143: }
144: VecRestoreArrayRead(x,&xx);
145: VecRestoreArrayRead(F,&FF);
146: return 0;
147: }
149: /*
151: Example function that when differenced produces the same matrix free Jacobian as FormFunction()
152: this is provided to show how a user can provide a different function
153: */
154: PetscErrorCode OtherFunctionForDifferencing(void *dummy,Vec x,Vec f)
155: {
158: FormFunction(NULL,x,f,dummy);
159: VecShift(f,1.0);
160: return 0;
161: }
163: /* -------------------- Form initial approximation ----------------- */
165: PetscErrorCode FormInitialGuess(SNES snes,Vec x)
166: {
168: PetscScalar pfive = .50;
169: VecSet(x,pfive);
170: return 0;
171: }
172: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
173: /* Evaluates a matrix that is used to precondition the matrix-free
174: jacobian. In this case, the explict preconditioner matrix is
175: also EXACTLY the Jacobian. In general, it would be some lower
176: order, simplified apprioximation */
178: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
179: {
180: const PetscScalar *xx;
181: PetscScalar A[3],d;
182: PetscInt i,n,j[3];
183: PetscErrorCode ierr;
184: AppCtx *user = (AppCtx*) dummy;
186: VecGetArrayRead(x,&xx);
187: VecGetSize(x,&n);
188: d = (PetscReal)(n - 1); d = d*d;
190: i = 0; A[0] = 1.0;
191: MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES);
192: for (i=1; i<n-1; i++) {
193: j[0] = i - 1; j[1] = i; j[2] = i + 1;
194: A[0] = d; A[1] = -2.0*d + 2.0*xx[i]; A[2] = d;
195: MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
196: }
197: i = n-1; A[0] = 1.0;
198: MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES);
199: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
200: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
201: VecRestoreArrayRead(x,&xx);
203: if (user->variant) {
204: MatMFFDSetBase(jac,x,NULL);
205: }
206: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
207: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
208: return 0;
209: }
211: PetscErrorCode FormJacobianNoMatrix(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
212: {
213: PetscErrorCode ierr;
214: AppCtx *user = (AppCtx*) dummy;
216: if (user->variant) {
217: MatMFFDSetBase(jac,x,NULL);
218: }
219: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
220: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
221: return 0;
222: }
224: /* -------------------- User-defined monitor ----------------------- */
226: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *dummy)
227: {
229: MonitorCtx *monP = (MonitorCtx*) dummy;
230: Vec x;
231: MPI_Comm comm;
233: PetscObjectGetComm((PetscObject)snes,&comm);
234: PetscFPrintf(comm,stdout,"iter = %D, SNES Function norm %g \n",its,(double)fnorm);
235: SNESGetSolution(snes,&x);
236: VecView(x,monP->viewer);
237: return 0;
238: }
241: /*TEST
243: test:
244: args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
246: test:
247: suffix: 2
248: args: -variant -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
249: output_file: output/ex7_1.out
251: # uses AIJ matrix to define diagonal matrix for Jacobian preconditioning
252: test:
253: suffix: 3
254: args: -variant -pc_type jacobi -snes_view -ksp_monitor
256: # uses MATMFFD matrix to define diagonal matrix for Jacobian preconditioning
257: test:
258: suffix: 4
259: args: -variant -pc_type jacobi -puremf -snes_view -ksp_monitor
261: TEST*/