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