Actual source code: ex24.c
petsc-3.3-p7 2013-05-11
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*,MatStructure*,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;
16: int main(int argc,char **argv)
17: {
18: TS ts; /* time integration context */
19: Vec X; /* solution, residual vectors */
20: Mat J; /* Jacobian matrix */
22: PetscScalar *x;
23: PetscReal ftime;
24: PetscInt i,steps,nits,lits;
25: PetscBool view_final;
26: Ctx ctx;
28: PetscInitialize(&argc,&argv,(char *)0,help);
29: ctx.n = 3;
30: PetscOptionsGetInt(PETSC_NULL,"-n",&ctx.n,PETSC_NULL);
31: if (ctx.n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The dimension specified with -n must be at least 2");
32: view_final = PETSC_FALSE;
33: PetscOptionsGetBool(PETSC_NULL,"-view_final",&view_final,PETSC_NULL);
34: ctx.monitor_short = PETSC_FALSE;
35: PetscOptionsGetBool(PETSC_NULL,"-monitor_short",&ctx.monitor_short,PETSC_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: MatGetVecs(J,&X,PETSC_NULL);
45: /* Create time integration context */
46: TSCreate(PETSC_COMM_WORLD,&ts);
47: TSSetType(ts,TSPSEUDO);
48: TSSetIFunction(ts,PETSC_NULL,FormIFunction,&ctx);
49: TSSetIJacobian(ts,J,J,FormIJacobian,&ctx);
50: TSSetDuration(ts,1000,1e14);
51: TSSetInitialTimeStep(ts,0.0,1e-3);
52: TSMonitorSet(ts,MonitorObjective,&ctx,PETSC_NULL);
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Customize time integrator; set runtime options
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: TSSetFromOptions(ts);
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Evaluate initial guess; then solve nonlinear system
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
62: VecSet(X,0.0);
63: VecGetArray(X,&x);
64: #if 1
65: x[0] = 5.;
66: x[1] = -5.;
67: for (i=2; i<ctx.n; i++) x[i] = 5.;
68: #else
69: x[0] = 1.0;
70: x[1] = 15.0;
71: for (i=2; i<ctx.n; i++) x[i] = 10.0;
72: #endif
73: VecRestoreArray(X,&x);
75: TSSolve(ts,X,&ftime);
76: TSGetTimeStepNumber(ts,&steps);
77: TSGetSNESIterations(ts,&nits);
78: TSGetKSPIterations(ts,&lits);
79: PetscPrintf(PETSC_COMM_WORLD,"Time integrator took (%D,%D,%D) iterations to reach final time %G\n",steps,nits,lits,ftime);
80: if (view_final) {
81: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
82: }
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Free work space. All PETSc objects should be destroyed when they
86: are no longer needed.
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: VecDestroy(&X);
90: MatDestroy(&J);
91: TSDestroy(&ts);
92: PetscFinalize();
93: return 0;
94: }
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++) {
114: f += PetscSqr(1. - x[i]) + 100. * PetscSqr(x[i+1] - PetscSqr(x[i]));
115: }
116: VecRestoreArrayRead(X,&x);
118: /* Compute norm of gradient */
119: VecDuplicate(X,&Xdot);
120: VecDuplicate(X,&F);
121: VecZeroEntries(Xdot);
122: FormIFunction(ts,t,X,Xdot,F,ictx);
123: VecNorm(F,NORM_2,&gnorm);
124: VecDestroy(&Xdot);
125: VecDestroy(&F);
127: TSGetTimeStep(ts,&dt);
128: TSGetSNES(ts,&snes);
129: SNESGetIterationNumber(snes,&snesit);
130: SNESGetLinearSolveIterations(snes,&linit);
131: PetscPrintf(PETSC_COMM_WORLD,
132: (ctx->monitor_short
133: ? "%3D t=%10.1e dt=%10.1e f=%10.1e df=%10.1e it=(%2D,%3D)\n"
134: : "%3D t=%10.4e dt=%10.4e f=%10.4e df=%10.4e it=(%2D,%3D)\n"),
135: step,(double)t,(double)dt,(double)PetscRealPart(f),(double)gnorm,snesit,linit);
136: return(0);
137: }
140: /* ------------------------------------------------------------------- */
143: /*
144: FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X))
146: Input Parameters:
147: + ts - the TS context
148: . t - time
149: . X - input vector
150: . Xdot - time derivative
151: - ctx - optional user-defined context
153: Output Parameter:
154: . F - function vector
155: */
156: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ictx)
157: {
159: const PetscScalar *x;
160: PetscScalar *f;
161: PetscInt i;
162: Ctx *ctx = (Ctx*)ictx;
165: /*
166: Get pointers to vector data.
167: - For default PETSc vectors, VecGetArray() returns a pointer to
168: the data array. Otherwise, the routine is implementation dependent.
169: - You MUST call VecRestoreArray() when you no longer need access to
170: the array.
171: */
172: VecGetArrayRead(X,&x);
173: VecZeroEntries(F);
174: VecGetArray(F,&f);
176: /* Compute gradient of objective */
177: for (i=0; i<ctx->n-1; i++) {
178: PetscScalar a,a0,a1;
179: a = x[i+1] - PetscSqr(x[i]);
180: a0 = -2.*x[i];
181: a1 = 1.;
182: f[i] += -2.*(1. - x[i]) + 200.*a*a0;
183: f[i+1] += 200.*a*a1;
184: }
185: /* Restore vectors */
186: VecRestoreArrayRead(X,&x);
187: VecRestoreArray(F,&f);
188: VecAXPY(F,1.0,Xdot);
189: return(0);
190: }
191: /* ------------------------------------------------------------------- */
194: /*
195: FormIJacobian - Evaluates Jacobian matrix.
197: Input Parameters:
198: + ts - the TS context
199: . t - pseudo-time
200: . X - input vector
201: . Xdot - time derivative
202: . shift - multiplier for mass matrix
203: . dummy - user-defined context
205: Output Parameters:
206: . J - Jacobian matrix
207: . B - optionally different preconditioning matrix
208: . flag - flag indicating matrix structure
209: */
210: static PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *J,Mat *B,MatStructure *flag,void *ictx)
211: {
212: const PetscScalar *x;
214: PetscInt i;
215: Ctx *ctx = (Ctx*)ictx;
218: MatZeroEntries(*B);
219: /*
220: Get pointer to vector data
221: */
222: VecGetArrayRead(X,&x);
224: /*
225: Compute Jacobian entries and insert into matrix.
226: */
227: for (i=0; i<ctx->n-1; i++) {
228: PetscInt rowcol[2];
229: PetscScalar v[2][2],a,a0,a1,a00,a01,a10,a11;
230: rowcol[0] = i;
231: rowcol[1] = i+1;
232: a = x[i+1] - PetscSqr(x[i]);
233: a0 = -2.*x[i];
234: a00 = -2.;
235: a01 = 0.;
236: a1 = 1.;
237: a10 = 0.;
238: a11 = 0.;
239: v[0][0] = 2. + 200.*(a*a00 + a0*a0);
240: v[0][1] = 200.*(a*a01 + a1*a0);
241: v[1][0] = 200.*(a*a10 + a0*a1);
242: v[1][1] = 200.*(a*a11 + a1*a1);
243: MatSetValues(*B,2,rowcol,2,rowcol,&v[0][0],ADD_VALUES);
244: }
245: for (i=0; i<ctx->n; i++) {
246: MatSetValue(*B,i,i,(PetscScalar)shift,ADD_VALUES);
247: }
249: VecRestoreArrayRead(X,&x);
251: /*
252: Assemble matrix
253: */
254: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
255: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
256: if (*J != *B){
257: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
258: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
259: }
260: return(0);
261: }