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