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