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 */
20: PetscScalar *x;
21: PetscReal ftime;
22: PetscInt i,steps,nits,lits;
23: PetscBool view_final;
24: Ctx ctx;
26: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
27: ctx.n = 3;
28: PetscOptionsGetInt(NULL,NULL,"-n",&ctx.n,NULL);
29: if (ctx.n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The dimension specified with -n must be at least 2");
31: view_final = PETSC_FALSE;
32: PetscOptionsGetBool(NULL,NULL,"-view_final",&view_final,NULL);
34: ctx.monitor_short = PETSC_FALSE;
35: PetscOptionsGetBool(NULL,NULL,"-monitor_short",&ctx.monitor_short,NULL);
37: /*
38: Create Jacobian matrix data structure and state vector
39: */
40: MatCreate(PETSC_COMM_WORLD,&J);
41: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,ctx.n,ctx.n);
42: MatSetFromOptions(J);
43: MatSetUp(J);
44: MatCreateVecs(J,&X,NULL);
46: /* Create time integration context */
47: TSCreate(PETSC_COMM_WORLD,&ts);
48: TSSetType(ts,TSPSEUDO);
49: TSSetIFunction(ts,NULL,FormIFunction,&ctx);
50: TSSetIJacobian(ts,J,J,FormIJacobian,&ctx);
51: TSSetMaxSteps(ts,1000);
52: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
53: TSSetTimeStep(ts,1e-3);
54: TSMonitorSet(ts,MonitorObjective,&ctx,NULL);
56: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: Customize time integrator; set runtime options
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: TSSetFromOptions(ts);
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Evaluate initial guess; then solve nonlinear system
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: VecSet(X,0.0);
65: VecGetArray(X,&x);
66: #if 1
67: x[0] = 5.;
68: x[1] = -5.;
69: for (i=2; i<ctx.n; i++) x[i] = 5.;
70: #else
71: x[0] = 1.0;
72: x[1] = 15.0;
73: for (i=2; i<ctx.n; i++) x[i] = 10.0;
74: #endif
75: VecRestoreArray(X,&x);
77: TSSolve(ts,X);
78: TSGetSolveTime(ts,&ftime);
79: TSGetStepNumber(ts,&steps);
80: TSGetSNESIterations(ts,&nits);
81: TSGetKSPIterations(ts,&lits);
82: PetscPrintf(PETSC_COMM_WORLD,"Time integrator took (%D,%D,%D) iterations to reach final time %g\n",steps,nits,lits,(double)ftime);
83: if (view_final) {
84: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
85: }
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Free work space. All PETSc objects should be destroyed when they
89: are no longer needed.
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: VecDestroy(&X);
93: MatDestroy(&J);
94: TSDestroy(&ts);
95: PetscFinalize();
96: return ierr;
97: }
99: static PetscErrorCode MonitorObjective(TS ts,PetscInt step,PetscReal t,Vec X,void *ictx)
100: {
101: Ctx *ctx = (Ctx*)ictx;
102: PetscErrorCode ierr;
103: const PetscScalar *x;
104: PetscScalar f;
105: PetscReal dt,gnorm;
106: PetscInt i,snesit,linit;
107: SNES snes;
108: Vec Xdot,F;
111: /* Compute objective functional */
112: VecGetArrayRead(X,&x);
113: f = 0;
114: for (i=0; i<ctx->n-1; i++) f += PetscSqr(1. - x[i]) + 100. * PetscSqr(x[i+1] - PetscSqr(x[i]));
115: VecRestoreArrayRead(X,&x);
117: /* Compute norm of gradient */
118: VecDuplicate(X,&Xdot);
119: VecDuplicate(X,&F);
120: VecZeroEntries(Xdot);
121: FormIFunction(ts,t,X,Xdot,F,ictx);
122: VecNorm(F,NORM_2,&gnorm);
123: VecDestroy(&Xdot);
124: VecDestroy(&F);
126: TSGetTimeStep(ts,&dt);
127: TSGetSNES(ts,&snes);
128: SNESGetIterationNumber(snes,&snesit);
129: SNESGetLinearSolveIterations(snes,&linit);
130: PetscPrintf(PETSC_COMM_WORLD,
131: (ctx->monitor_short
132: ? "%3D t=%10.1e dt=%10.1e f=%10.1e df=%10.1e it=(%2D,%3D)\n"
133: : "%3D t=%10.4e dt=%10.4e f=%10.4e df=%10.4e it=(%2D,%3D)\n"),
134: step,(double)t,(double)dt,(double)PetscRealPart(f),(double)gnorm,snesit,linit);
135: return(0);
136: }
138: /* ------------------------------------------------------------------- */
139: /*
140: FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X))
142: Input Parameters:
143: + ts - the TS context
144: . t - time
145: . X - input vector
146: . Xdot - time derivative
147: - ctx - optional user-defined context
149: Output Parameter:
150: . F - function vector
151: */
152: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ictx)
153: {
154: PetscErrorCode ierr;
155: const PetscScalar *x;
156: PetscScalar *f;
157: PetscInt i;
158: Ctx *ctx = (Ctx*)ictx;
161: /*
162: Get pointers to vector data.
163: - For default PETSc vectors, VecGetArray() returns a pointer to
164: the data array. Otherwise, the routine is implementation dependent.
165: - You MUST call VecRestoreArray() when you no longer need access to
166: the array.
167: */
168: VecGetArrayRead(X,&x);
169: VecZeroEntries(F);
170: VecGetArray(F,&f);
172: /* Compute gradient of objective */
173: for (i=0; i<ctx->n-1; i++) {
174: PetscScalar a,a0,a1;
175: a = x[i+1] - PetscSqr(x[i]);
176: a0 = -2.*x[i];
177: a1 = 1.;
178: f[i] += -2.*(1. - x[i]) + 200.*a*a0;
179: f[i+1] += 200.*a*a1;
180: }
181: /* Restore vectors */
182: VecRestoreArrayRead(X,&x);
183: VecRestoreArray(F,&f);
184: VecAXPY(F,1.0,Xdot);
185: return(0);
186: }
187: /* ------------------------------------------------------------------- */
188: /*
189: FormIJacobian - Evaluates Jacobian matrix.
191: Input Parameters:
192: + ts - the TS context
193: . t - pseudo-time
194: . X - input vector
195: . Xdot - time derivative
196: . shift - multiplier for mass matrix
197: . dummy - user-defined context
199: Output Parameters:
200: . J - Jacobian matrix
201: . B - optionally different preconditioning matrix
202: . flag - flag indicating matrix structure
203: */
204: static PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat J,Mat B,void *ictx)
205: {
206: const PetscScalar *x;
207: PetscErrorCode ierr;
208: PetscInt i;
209: Ctx *ctx = (Ctx*)ictx;
212: MatZeroEntries(B);
213: /*
214: Get pointer to vector data
215: */
216: VecGetArrayRead(X,&x);
218: /*
219: Compute Jacobian entries and insert into matrix.
220: */
221: for (i=0; i<ctx->n-1; i++) {
222: PetscInt rowcol[2];
223: PetscScalar v[2][2],a,a0,a1,a00,a01,a10,a11;
224: rowcol[0] = i;
225: rowcol[1] = i+1;
226: a = x[i+1] - PetscSqr(x[i]);
227: a0 = -2.*x[i];
228: a00 = -2.;
229: a01 = 0.;
230: a1 = 1.;
231: a10 = 0.;
232: a11 = 0.;
233: v[0][0] = 2. + 200.*(a*a00 + a0*a0);
234: v[0][1] = 200.*(a*a01 + a1*a0);
235: v[1][0] = 200.*(a*a10 + a0*a1);
236: v[1][1] = 200.*(a*a11 + a1*a1);
237: MatSetValues(B,2,rowcol,2,rowcol,&v[0][0],ADD_VALUES);
238: }
239: for (i=0; i<ctx->n; i++) {
240: MatSetValue(B,i,i,(PetscScalar)shift,ADD_VALUES);
241: }
243: VecRestoreArrayRead(X,&x);
245: /*
246: Assemble matrix
247: */
248: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
249: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
250: if (J != B) {
251: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
252: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
253: }
254: return(0);
255: }
257: /*TEST
259: test:
260: requires: !single
262: test:
263: 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
264: requires: !single
265: suffix: 2
267: TEST*/