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, NULL, 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" : "%3" PetscInt_FMT " t=%10.4e dt=%10.4e f=%10.4e df=%10.4e it=(%2" PetscInt_FMT ",%3" PetscInt_FMT ")\n", step, (double)t, (double)dt, (double)PetscRealPart(f), (double)gnorm, snesit, linit)); 128 PetscFunctionReturn(PETSC_SUCCESS); 129 } 130 131 /* ------------------------------------------------------------------- */ 132 /* 133 FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X)) 134 135 Input Parameters: 136 + ts - the TS context 137 . t - time 138 . X - input vector 139 . Xdot - time derivative 140 - ctx - optional user-defined context 141 142 Output Parameter: 143 . F - function vector 144 */ 145 static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ictx) 146 { 147 const PetscScalar *x; 148 PetscScalar *f; 149 PetscInt i; 150 Ctx *ctx = (Ctx *)ictx; 151 152 PetscFunctionBeginUser; 153 /* 154 Get pointers to vector data. 155 - For default PETSc vectors, VecGetArray() returns a pointer to 156 the data array. Otherwise, the routine is implementation dependent. 157 - You MUST call VecRestoreArray() when you no longer need access to 158 the array. 159 */ 160 PetscCall(VecGetArrayRead(X, &x)); 161 PetscCall(VecZeroEntries(F)); 162 PetscCall(VecGetArray(F, &f)); 163 164 /* Compute gradient of objective */ 165 for (i = 0; i < ctx->n - 1; i++) { 166 PetscScalar a, a0, a1; 167 a = x[i + 1] - PetscSqr(x[i]); 168 a0 = -2. * x[i]; 169 a1 = 1.; 170 f[i] += -2. * (1. - x[i]) + 200. * a * a0; 171 f[i + 1] += 200. * a * a1; 172 } 173 /* Restore vectors */ 174 PetscCall(VecRestoreArrayRead(X, &x)); 175 PetscCall(VecRestoreArray(F, &f)); 176 PetscCall(VecAXPY(F, 1.0, Xdot)); 177 PetscFunctionReturn(PETSC_SUCCESS); 178 } 179 /* ------------------------------------------------------------------- */ 180 /* 181 FormIJacobian - Evaluates Jacobian matrix. 182 183 Input Parameters: 184 + ts - the TS context 185 . t - pseudo-time 186 . X - input vector 187 . Xdot - time derivative 188 . shift - multiplier for mass matrix 189 . dummy - user-defined context 190 191 Output Parameters: 192 . J - Jacobian matrix 193 . B - optionally different preconditioning matrix 194 */ 195 static PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat J, Mat B, void *ictx) 196 { 197 const PetscScalar *x; 198 PetscInt i; 199 Ctx *ctx = (Ctx *)ictx; 200 201 PetscFunctionBeginUser; 202 PetscCall(MatZeroEntries(B)); 203 /* 204 Get pointer to vector data 205 */ 206 PetscCall(VecGetArrayRead(X, &x)); 207 208 /* 209 Compute Jacobian entries and insert into matrix. 210 */ 211 for (i = 0; i < ctx->n - 1; i++) { 212 PetscInt rowcol[2]; 213 PetscScalar v[2][2], a, a0, a1, a00, a01, a10, a11; 214 rowcol[0] = i; 215 rowcol[1] = i + 1; 216 a = x[i + 1] - PetscSqr(x[i]); 217 a0 = -2. * x[i]; 218 a00 = -2.; 219 a01 = 0.; 220 a1 = 1.; 221 a10 = 0.; 222 a11 = 0.; 223 v[0][0] = 2. + 200. * (a * a00 + a0 * a0); 224 v[0][1] = 200. * (a * a01 + a1 * a0); 225 v[1][0] = 200. * (a * a10 + a0 * a1); 226 v[1][1] = 200. * (a * a11 + a1 * a1); 227 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &v[0][0], ADD_VALUES)); 228 } 229 for (i = 0; i < ctx->n; i++) PetscCall(MatSetValue(B, i, i, (PetscScalar)shift, ADD_VALUES)); 230 231 PetscCall(VecRestoreArrayRead(X, &x)); 232 233 /* 234 Assemble matrix 235 */ 236 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 237 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 238 if (J != B) { 239 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 240 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 241 } 242 PetscFunctionReturn(PETSC_SUCCESS); 243 } 244 245 /*TEST 246 247 test: 248 requires: !single 249 250 test: 251 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 unlimited 252 requires: !single 253 suffix: 2 254 255 TEST*/ 256