1 static char help[] = "Time-dependent PDE in 2d for calculating joint PDF. \n"; 2 /* 3 p_t = -x_t*p_x -y_t*p_y + f(t)*p_yy 4 xmin < x < xmax, ymin < y < ymax; 5 6 Boundary conditions: 7 Zero dirichlet in y using ghosted values 8 Periodic in x 9 10 Note that x_t and y_t in the above are given functions of x and y; they are not derivatives of x and y. 11 x_t = (y - ws) 12 y_t = (ws/2H)*(Pm - Pmax*sin(x) - D*(w - ws)) 13 14 In this example, we can see the effect of a fault, that zeroes the electrical power output 15 Pmax*sin(x), on the PDF. The fault on/off times can be controlled by options -tf and -tcl respectively. 16 17 */ 18 19 #include <petscdm.h> 20 #include <petscdmda.h> 21 #include <petscts.h> 22 23 /* 24 User-defined data structures and routines 25 */ 26 typedef struct { 27 PetscScalar ws; /* Synchronous speed */ 28 PetscScalar H; /* Inertia constant */ 29 PetscScalar D; /* Damping constant */ 30 PetscScalar Pmax, Pmax_s; /* Maximum power output of generator */ 31 PetscScalar PM_min; /* Mean mechanical power input */ 32 PetscScalar lambda; /* correlation time */ 33 PetscScalar q; /* noise strength */ 34 PetscScalar mux; /* Initial average angle */ 35 PetscScalar sigmax; /* Standard deviation of initial angle */ 36 PetscScalar muy; /* Average speed */ 37 PetscScalar sigmay; /* standard deviation of initial speed */ 38 PetscScalar rho; /* Cross-correlation coefficient */ 39 PetscScalar xmin; /* left boundary of angle */ 40 PetscScalar xmax; /* right boundary of angle */ 41 PetscScalar ymin; /* bottom boundary of speed */ 42 PetscScalar ymax; /* top boundary of speed */ 43 PetscScalar dx; /* x step size */ 44 PetscScalar dy; /* y step size */ 45 PetscScalar disper_coe; /* Dispersion coefficient */ 46 DM da; 47 PetscInt st_width; /* Stencil width */ 48 DMBoundaryType bx; /* x boundary type */ 49 DMBoundaryType by; /* y boundary type */ 50 PetscReal tf, tcl; /* Fault incidence and clearing times */ 51 } AppCtx; 52 53 PetscErrorCode Parameter_settings(AppCtx *); 54 PetscErrorCode ini_bou(Vec, AppCtx *); 55 PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *); 56 PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *); 57 PetscErrorCode PostStep(TS); 58 59 int main(int argc, char **argv) 60 { 61 Vec x; /* Solution vector */ 62 TS ts; /* Time-stepping context */ 63 AppCtx user; /* Application context */ 64 PetscViewer viewer; 65 66 PetscFunctionBeginUser; 67 PetscCall(PetscInitialize(&argc, &argv, "petscopt_ex8", help)); 68 69 /* Get physics and time parameters */ 70 PetscCall(Parameter_settings(&user)); 71 /* Create a 2D DA with dof = 1 */ 72 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, user.bx, user.by, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, user.st_width, NULL, NULL, &user.da)); 73 PetscCall(DMSetFromOptions(user.da)); 74 PetscCall(DMSetUp(user.da)); 75 /* Set x and y coordinates */ 76 PetscCall(DMDASetUniformCoordinates(user.da, user.xmin, user.xmax, user.ymin, user.ymax, 0, 0)); 77 PetscCall(DMDASetCoordinateName(user.da, 0, "X - the angle")); 78 PetscCall(DMDASetCoordinateName(user.da, 1, "Y - the speed")); 79 80 /* Get global vector x from DM */ 81 PetscCall(DMCreateGlobalVector(user.da, &x)); 82 83 PetscCall(ini_bou(x, &user)); 84 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "ini_x", FILE_MODE_WRITE, &viewer)); 85 PetscCall(VecView(x, viewer)); 86 PetscCall(PetscViewerDestroy(&viewer)); 87 88 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 89 PetscCall(TSSetDM(ts, user.da)); 90 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 91 PetscCall(TSSetType(ts, TSARKIMEX)); 92 PetscCall(TSSetIFunction(ts, NULL, IFunction, &user)); 93 /* PetscCall(TSSetIJacobian(ts,NULL,NULL,IJacobian,&user)); */ 94 PetscCall(TSSetApplicationContext(ts, &user)); 95 PetscCall(TSSetTimeStep(ts, .005)); 96 PetscCall(TSSetFromOptions(ts)); 97 PetscCall(TSSetPostStep(ts, PostStep)); 98 PetscCall(TSSolve(ts, x)); 99 100 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "fin_x", FILE_MODE_WRITE, &viewer)); 101 PetscCall(VecView(x, viewer)); 102 PetscCall(PetscViewerDestroy(&viewer)); 103 104 PetscCall(VecDestroy(&x)); 105 PetscCall(DMDestroy(&user.da)); 106 PetscCall(TSDestroy(&ts)); 107 PetscCall(PetscFinalize()); 108 return 0; 109 } 110 111 PetscErrorCode PostStep(TS ts) 112 { 113 Vec X; 114 AppCtx *user; 115 PetscReal t; 116 PetscScalar asum; 117 118 PetscFunctionBegin; 119 PetscCall(TSGetApplicationContext(ts, &user)); 120 PetscCall(TSGetTime(ts, &t)); 121 PetscCall(TSGetSolution(ts, &X)); 122 /* 123 if (t >= .2) { 124 PetscCall(TSGetSolution(ts,&X)); 125 PetscCall(VecView(X,PETSC_VIEWER_BINARY_WORLD)); 126 exit(0); 127 results in initial conditions after fault in binaryoutput 128 }*/ 129 130 if ((t > user->tf) && (t < user->tcl)) user->Pmax = 0.0; /* A short-circuit that drives the electrical power output (Pmax*sin(delta)) to zero */ 131 else user->Pmax = user->Pmax_s; 132 133 PetscCall(VecSum(X, &asum)); 134 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sum(p) at t = %f = %f\n", (double)t, (double)(asum))); 135 PetscFunctionReturn(PETSC_SUCCESS); 136 } 137 138 PetscErrorCode ini_bou(Vec X, AppCtx *user) 139 { 140 DM cda; 141 DMDACoor2d **coors; 142 PetscScalar **p; 143 Vec gc; 144 PetscInt M, N, Ir, J; 145 PetscMPIInt rank; 146 147 PetscFunctionBeginUser; 148 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 149 PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 150 user->dx = (user->xmax - user->xmin) / (M - 1); 151 user->dy = (user->ymax - user->ymin) / (N - 1); 152 PetscCall(DMGetCoordinateDM(user->da, &cda)); 153 PetscCall(DMGetCoordinates(user->da, &gc)); 154 PetscCall(DMDAVecGetArrayRead(cda, gc, &coors)); 155 PetscCall(DMDAVecGetArray(user->da, X, &p)); 156 157 /* Point mass at (mux,muy) */ 158 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original user->mux = %f, user->muy = %f\n", (double)PetscRealPart(user->mux), (double)PetscRealPart(user->muy))); 159 PetscCall(DMDAGetLogicalCoordinate(user->da, user->mux, user->muy, 0.0, &Ir, &J, NULL, &user->mux, &user->muy, NULL)); 160 user->PM_min = user->Pmax * PetscSinScalar(user->mux); 161 PetscCall( 162 PetscPrintf(PETSC_COMM_WORLD, "Corrected user->mux = %f, user->muy = %f user->PM_min = %f,user->dx = %f\n", (double)PetscRealPart(user->mux), (double)PetscRealPart(user->muy), (double)PetscRealPart(user->PM_min), (double)PetscRealPart(user->dx))); 163 if (Ir > -1 && J > -1) p[J][Ir] = 1.0; 164 165 PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors)); 166 PetscCall(DMDAVecRestoreArray(user->da, X, &p)); 167 PetscFunctionReturn(PETSC_SUCCESS); 168 } 169 170 /* First advection term */ 171 PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar *p1, AppCtx *user) 172 { 173 PetscScalar f, fpos, fneg; 174 PetscFunctionBegin; 175 f = (y - user->ws); 176 fpos = PetscMax(f, 0); 177 fneg = PetscMin(f, 0); 178 if (user->st_width == 1) { 179 *p1 = fpos * (p[j][i] - p[j][i - 1]) / user->dx + fneg * (p[j][i + 1] - p[j][i]) / user->dx; 180 } else if (user->st_width == 2) { 181 *p1 = fpos * (3 * p[j][i] - 4 * p[j][i - 1] + p[j][i - 2]) / (2 * user->dx) + fneg * (-p[j][i + 2] + 4 * p[j][i + 1] - 3 * p[j][i]) / (2 * user->dx); 182 } else if (user->st_width == 3) { 183 *p1 = fpos * (2 * p[j][i + 1] + 3 * p[j][i] - 6 * p[j][i - 1] + p[j][i - 2]) / (6 * user->dx) + fneg * (-p[j][i + 2] + 6 * p[j][i + 1] - 3 * p[j][i] - 2 * p[j][i - 1]) / (6 * user->dx); 184 } 185 /* *p1 = f*(p[j][i+1] - p[j][i-1])/user->dx;*/ 186 PetscFunctionReturn(PETSC_SUCCESS); 187 } 188 189 /* Second advection term */ 190 PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscScalar y, PetscInt i, PetscInt j, PetscInt N, PetscScalar *p2, AppCtx *user) 191 { 192 PetscScalar f, fpos, fneg; 193 PetscFunctionBegin; 194 f = (user->ws / (2 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(x) - user->D * (y - user->ws)); 195 fpos = PetscMax(f, 0); 196 fneg = PetscMin(f, 0); 197 if (user->st_width == 1) { 198 *p2 = fpos * (p[j][i] - p[j - 1][i]) / user->dy + fneg * (p[j + 1][i] - p[j][i]) / user->dy; 199 } else if (user->st_width == 2) { 200 *p2 = fpos * (3 * p[j][i] - 4 * p[j - 1][i] + p[j - 2][i]) / (2 * user->dy) + fneg * (-p[j + 2][i] + 4 * p[j + 1][i] - 3 * p[j][i]) / (2 * user->dy); 201 } else if (user->st_width == 3) { 202 *p2 = fpos * (2 * p[j + 1][i] + 3 * p[j][i] - 6 * p[j - 1][i] + p[j - 2][i]) / (6 * user->dy) + fneg * (-p[j + 2][i] + 6 * p[j + 1][i] - 3 * p[j][i] - 2 * p[j - 1][i]) / (6 * user->dy); 203 } 204 205 /* *p2 = f*(p[j+1][i] - p[j-1][i])/user->dy;*/ 206 PetscFunctionReturn(PETSC_SUCCESS); 207 } 208 209 /* Diffusion term */ 210 PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user) 211 { 212 PetscFunctionBeginUser; 213 if (user->st_width == 1) { 214 *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy)); 215 } else if (user->st_width == 2) { 216 *p_diff = user->disper_coe * ((-p[j - 2][i] + 16 * p[j - 1][i] - 30 * p[j][i] + 16 * p[j + 1][i] - p[j + 2][i]) / (12.0 * user->dy * user->dy)); 217 } else if (user->st_width == 3) { 218 *p_diff = user->disper_coe * ((2 * p[j - 3][i] - 27 * p[j - 2][i] + 270 * p[j - 1][i] - 490 * p[j][i] + 270 * p[j + 1][i] - 27 * p[j + 2][i] + 2 * p[j + 3][i]) / (180.0 * user->dy * user->dy)); 219 } 220 PetscFunctionReturn(PETSC_SUCCESS); 221 } 222 223 PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx) 224 { 225 AppCtx *user = (AppCtx *)ctx; 226 DM cda; 227 DMDACoor2d **coors; 228 PetscScalar **p, **f, **pdot; 229 PetscInt i, j; 230 PetscInt xs, ys, xm, ym, M, N; 231 Vec localX, gc, localXdot; 232 PetscScalar p_adv1 = 0.0, p_adv2 = 0.0, p_diff = 0; 233 PetscScalar diffuse1, gamma; 234 235 PetscFunctionBeginUser; 236 PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 237 PetscCall(DMGetCoordinateDM(user->da, &cda)); 238 PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0)); 239 240 PetscCall(DMGetLocalVector(user->da, &localX)); 241 PetscCall(DMGetLocalVector(user->da, &localXdot)); 242 243 PetscCall(DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX)); 244 PetscCall(DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX)); 245 PetscCall(DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot)); 246 PetscCall(DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot)); 247 248 PetscCall(DMGetCoordinatesLocal(user->da, &gc)); 249 250 PetscCall(DMDAVecGetArrayRead(cda, gc, &coors)); 251 PetscCall(DMDAVecGetArrayRead(user->da, localX, &p)); 252 PetscCall(DMDAVecGetArrayRead(user->da, localXdot, &pdot)); 253 PetscCall(DMDAVecGetArray(user->da, F, &f)); 254 255 gamma = user->D * user->ws / (2 * user->H); 256 diffuse1 = user->lambda * user->lambda * user->q / (user->lambda * gamma + 1) * (1.0 - PetscExpScalar(-t * (gamma + 1.0) / user->lambda)); 257 user->disper_coe = user->ws * user->ws / (4 * user->H * user->H) * diffuse1; 258 259 for (i = xs; i < xs + xm; i++) { 260 for (j = ys; j < ys + ym; j++) { 261 PetscCall(adv1(p, coors[j][i].y, i, j, M, &p_adv1, user)); 262 PetscCall(adv2(p, coors[j][i].x, coors[j][i].y, i, j, N, &p_adv2, user)); 263 PetscCall(diffuse(p, i, j, t, &p_diff, user)); 264 f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i]; 265 } 266 } 267 PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &p)); 268 PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &pdot)); 269 PetscCall(DMRestoreLocalVector(user->da, &localX)); 270 PetscCall(DMRestoreLocalVector(user->da, &localXdot)); 271 PetscCall(DMDAVecRestoreArray(user->da, F, &f)); 272 PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors)); 273 274 PetscFunctionReturn(PETSC_SUCCESS); 275 } 276 277 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ctx) 278 { 279 AppCtx *user = (AppCtx *)ctx; 280 DM cda; 281 DMDACoor2d **coors; 282 PetscInt i, j; 283 PetscInt xs, ys, xm, ym, M, N; 284 Vec gc; 285 PetscScalar val[5], xi, yi; 286 MatStencil row, col[5]; 287 PetscScalar c1, c3, c5, c1pos, c1neg, c3pos, c3neg; 288 289 PetscFunctionBeginUser; 290 PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 291 PetscCall(DMGetCoordinateDM(user->da, &cda)); 292 PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0)); 293 294 PetscCall(DMGetCoordinatesLocal(user->da, &gc)); 295 PetscCall(DMDAVecGetArrayRead(cda, gc, &coors)); 296 for (i = xs; i < xs + xm; i++) { 297 for (j = ys; j < ys + ym; j++) { 298 PetscInt nc = 0; 299 xi = coors[j][i].x; 300 yi = coors[j][i].y; 301 row.i = i; 302 row.j = j; 303 c1 = (yi - user->ws) / user->dx; 304 c1pos = PetscMax(c1, 0); 305 c1neg = PetscMin(c1, 0); 306 c3 = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi) - user->D * (yi - user->ws)) / user->dy; 307 c3pos = PetscMax(c3, 0); 308 c3neg = PetscMin(c3, 0); 309 c5 = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy); 310 col[nc].i = i - 1; 311 col[nc].j = j; 312 val[nc++] = c1pos; 313 col[nc].i = i + 1; 314 col[nc].j = j; 315 val[nc++] = -c1neg; 316 col[nc].i = i; 317 col[nc].j = j - 1; 318 val[nc++] = c3pos + c5; 319 col[nc].i = i; 320 col[nc].j = j + 1; 321 val[nc++] = -c3neg + c5; 322 col[nc].i = i; 323 col[nc].j = j; 324 val[nc++] = -c1pos + c1neg - c3pos + c3neg - 2 * c5 - a; 325 PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES)); 326 } 327 } 328 PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors)); 329 330 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY)); 331 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY)); 332 if (J != Jpre) { 333 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 334 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 335 } 336 PetscFunctionReturn(PETSC_SUCCESS); 337 } 338 339 PetscErrorCode Parameter_settings(AppCtx *user) 340 { 341 PetscBool flg; 342 343 PetscFunctionBeginUser; 344 /* Set default parameters */ 345 user->ws = 1.0; 346 user->H = 5.0; 347 user->D = 0.0; 348 user->Pmax = user->Pmax_s = 2.1; 349 user->PM_min = 1.0; 350 user->lambda = 0.1; 351 user->q = 1.0; 352 user->mux = PetscAsinScalar(user->PM_min / user->Pmax); 353 user->sigmax = 0.1; 354 user->sigmay = 0.1; 355 user->rho = 0.0; 356 user->xmin = -PETSC_PI; 357 user->xmax = PETSC_PI; 358 user->bx = DM_BOUNDARY_PERIODIC; 359 user->by = DM_BOUNDARY_GHOSTED; 360 user->tf = user->tcl = -1; 361 user->ymin = -2.0; 362 user->ymax = 2.0; 363 user->st_width = 1; 364 365 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg)); 366 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg)); 367 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-D", &user->D, &flg)); 368 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg)); 369 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg)); 370 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg)); 371 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg)); 372 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg)); 373 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg)); 374 if (flg == 0) user->muy = user->ws; 375 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg)); 376 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg)); 377 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg)); 378 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg)); 379 PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_width", &user->st_width, &flg)); 380 PetscCall(PetscOptionsGetEnum(NULL, NULL, "-bx", DMBoundaryTypes, (PetscEnum *)&user->bx, &flg)); 381 PetscCall(PetscOptionsGetEnum(NULL, NULL, "-by", DMBoundaryTypes, (PetscEnum *)&user->by, &flg)); 382 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tf", &user->tf, &flg)); 383 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tcl", &user->tcl, &flg)); 384 PetscFunctionReturn(PETSC_SUCCESS); 385 } 386 387 /*TEST 388 389 build: 390 requires: !complex x 391 392 test: 393 args: -ts_max_steps 1 394 localrunfiles: petscopt_ex8 395 timeoutfactor: 3 396 397 TEST*/ 398