Actual source code: ex24.c
1: static char help[] = "Pseudotransient continuation to solve a many-variable system that comes from the 2 variable Rosenbrock function + trivial.\n\n";
3: #include <petscts.h>
5: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
6: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
7: static PetscErrorCode MonitorObjective(TS,PetscInt,PetscReal,Vec,void*);
9: typedef struct {
10: PetscInt n;
11: PetscBool monitor_short;
12: } Ctx;
14: int main(int argc,char **argv)
15: {
16: TS ts; /* time integration context */
17: Vec X; /* solution, residual vectors */
18: Mat J; /* Jacobian matrix */
19: PetscScalar *x;
20: PetscReal ftime;
21: PetscInt i,steps,nits,lits;
22: PetscBool view_final;
23: Ctx ctx;
25: PetscInitialize(&argc,&argv,(char*)0,help);
26: ctx.n = 3;
27: PetscOptionsGetInt(NULL,NULL,"-n",&ctx.n,NULL);
30: view_final = PETSC_FALSE;
31: PetscOptionsGetBool(NULL,NULL,"-view_final",&view_final,NULL);
33: ctx.monitor_short = PETSC_FALSE;
34: PetscOptionsGetBool(NULL,NULL,"-monitor_short",&ctx.monitor_short,NULL);
36: /*
37: Create Jacobian matrix data structure and state vector
38: */
39: MatCreate(PETSC_COMM_WORLD,&J);
40: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,ctx.n,ctx.n);
41: MatSetFromOptions(J);
42: MatSetUp(J);
43: MatCreateVecs(J,&X,NULL);
45: /* Create time integration context */
46: TSCreate(PETSC_COMM_WORLD,&ts);
47: TSSetType(ts,TSPSEUDO);
48: TSSetIFunction(ts,NULL,FormIFunction,&ctx);
49: TSSetIJacobian(ts,J,J,FormIJacobian,&ctx);
50: TSSetMaxSteps(ts,1000);
51: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
52: TSSetTimeStep(ts,1e-3);
53: TSMonitorSet(ts,MonitorObjective,&ctx,NULL);
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Customize time integrator; set runtime options
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58: TSSetFromOptions(ts);
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Evaluate initial guess; then solve nonlinear system
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: VecSet(X,0.0);
64: VecGetArray(X,&x);
65: #if 1
66: x[0] = 5.;
67: x[1] = -5.;
68: for (i=2; i<ctx.n; i++) x[i] = 5.;
69: #else
70: x[0] = 1.0;
71: x[1] = 15.0;
72: for (i=2; i<ctx.n; i++) x[i] = 10.0;
73: #endif
74: VecRestoreArray(X,&x);
76: TSSolve(ts,X);
77: TSGetSolveTime(ts,&ftime);
78: TSGetStepNumber(ts,&steps);
79: TSGetSNESIterations(ts,&nits);
80: TSGetKSPIterations(ts,&lits);
81: PetscPrintf(PETSC_COMM_WORLD,"Time integrator took (%D,%D,%D) iterations to reach final time %g\n",steps,nits,lits,(double)ftime);
82: if (view_final) {
83: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
84: }
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Free work space. All PETSc objects should be destroyed when they
88: are no longer needed.
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: VecDestroy(&X);
92: MatDestroy(&J);
93: TSDestroy(&ts);
94: PetscFinalize();
95: return 0;
96: }
98: static PetscErrorCode MonitorObjective(TS ts,PetscInt step,PetscReal t,Vec X,void *ictx)
99: {
100: Ctx *ctx = (Ctx*)ictx;
101: PetscErrorCode ierr;
102: const PetscScalar *x;
103: PetscScalar f;
104: PetscReal dt,gnorm;
105: PetscInt i,snesit,linit;
106: SNES snes;
107: Vec Xdot,F;
110: /* Compute objective functional */
111: VecGetArrayRead(X,&x);
112: f = 0;
113: for (i=0; i<ctx->n-1; i++) f += PetscSqr(1. - x[i]) + 100. * PetscSqr(x[i+1] - PetscSqr(x[i]));
114: VecRestoreArrayRead(X,&x);
116: /* Compute norm of gradient */
117: VecDuplicate(X,&Xdot);
118: VecDuplicate(X,&F);
119: VecZeroEntries(Xdot);
120: FormIFunction(ts,t,X,Xdot,F,ictx);
121: VecNorm(F,NORM_2,&gnorm);
122: VecDestroy(&Xdot);
123: VecDestroy(&F);
125: TSGetTimeStep(ts,&dt);
126: TSGetSNES(ts,&snes);
127: SNESGetIterationNumber(snes,&snesit);
128: SNESGetLinearSolveIterations(snes,&linit);
129: PetscPrintf(PETSC_COMM_WORLD,
130: (ctx->monitor_short
131: ? "%3D t=%10.1e dt=%10.1e f=%10.1e df=%10.1e it=(%2D,%3D)\n"
132: : "%3D t=%10.4e dt=%10.4e f=%10.4e df=%10.4e it=(%2D,%3D)\n"),
133: step,(double)t,(double)dt,(double)PetscRealPart(f),(double)gnorm,snesit,linit);
134: return 0;
135: }
137: /* ------------------------------------------------------------------- */
138: /*
139: FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X))
141: Input Parameters:
142: + ts - the TS context
143: . t - time
144: . X - input vector
145: . Xdot - time derivative
146: - ctx - optional user-defined context
148: Output Parameter:
149: . F - function vector
150: */
151: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ictx)
152: {
153: const PetscScalar *x;
154: PetscScalar *f;
155: PetscInt i;
156: Ctx *ctx = (Ctx*)ictx;
159: /*
160: Get pointers to vector data.
161: - For default PETSc vectors, VecGetArray() returns a pointer to
162: the data array. Otherwise, the routine is implementation dependent.
163: - You MUST call VecRestoreArray() when you no longer need access to
164: the array.
165: */
166: VecGetArrayRead(X,&x);
167: VecZeroEntries(F);
168: VecGetArray(F,&f);
170: /* Compute gradient of objective */
171: for (i=0; i<ctx->n-1; i++) {
172: PetscScalar a,a0,a1;
173: a = x[i+1] - PetscSqr(x[i]);
174: a0 = -2.*x[i];
175: a1 = 1.;
176: f[i] += -2.*(1. - x[i]) + 200.*a*a0;
177: f[i+1] += 200.*a*a1;
178: }
179: /* Restore vectors */
180: VecRestoreArrayRead(X,&x);
181: VecRestoreArray(F,&f);
182: VecAXPY(F,1.0,Xdot);
183: return 0;
184: }
185: /* ------------------------------------------------------------------- */
186: /*
187: FormIJacobian - Evaluates Jacobian matrix.
189: Input Parameters:
190: + ts - the TS context
191: . t - pseudo-time
192: . X - input vector
193: . Xdot - time derivative
194: . shift - multiplier for mass matrix
195: . dummy - user-defined context
197: Output Parameters:
198: . J - Jacobian matrix
199: . B - optionally different preconditioning matrix
200: . flag - flag indicating matrix structure
201: */
202: static PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat J,Mat B,void *ictx)
203: {
204: const PetscScalar *x;
205: PetscInt i;
206: Ctx *ctx = (Ctx*)ictx;
209: MatZeroEntries(B);
210: /*
211: Get pointer to vector data
212: */
213: VecGetArrayRead(X,&x);
215: /*
216: Compute Jacobian entries and insert into matrix.
217: */
218: for (i=0; i<ctx->n-1; i++) {
219: PetscInt rowcol[2];
220: PetscScalar v[2][2],a,a0,a1,a00,a01,a10,a11;
221: rowcol[0] = i;
222: rowcol[1] = i+1;
223: a = x[i+1] - PetscSqr(x[i]);
224: a0 = -2.*x[i];
225: a00 = -2.;
226: a01 = 0.;
227: a1 = 1.;
228: a10 = 0.;
229: a11 = 0.;
230: v[0][0] = 2. + 200.*(a*a00 + a0*a0);
231: v[0][1] = 200.*(a*a01 + a1*a0);
232: v[1][0] = 200.*(a*a10 + a0*a1);
233: v[1][1] = 200.*(a*a11 + a1*a1);
234: MatSetValues(B,2,rowcol,2,rowcol,&v[0][0],ADD_VALUES);
235: }
236: for (i=0; i<ctx->n; i++) {
237: MatSetValue(B,i,i,(PetscScalar)shift,ADD_VALUES);
238: }
240: VecRestoreArrayRead(X,&x);
242: /*
243: Assemble matrix
244: */
245: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
246: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
247: if (J != B) {
248: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
249: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
250: }
251: return 0;
252: }
254: /*TEST
256: test:
257: requires: !single
259: test:
260: args: -pc_type lu -ts_dt 1e-5 -ts_max_time 1e5 -n 50 -monitor_short -snes_max_it 5 -snes_type newtonls -ts_max_snes_failures -1
261: requires: !single
262: suffix: 2
264: TEST*/