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 x_t = (y - ws) y_t = (ws/2H)*(Pm - Pmax*sin(x)) 6 7 Boundary conditions: -bc_type 0 => Zero dirichlet boundary 8 -bc_type 1 => Steady state boundary condition 9 Steady state boundary condition found by setting p_t = 0 10 */ 11 12 #include <petscdm.h> 13 #include <petscdmda.h> 14 #include <petscts.h> 15 16 /* 17 User-defined data structures and routines 18 */ 19 typedef struct { 20 PetscScalar ws; /* Synchronous speed */ 21 PetscScalar H; /* Inertia constant */ 22 PetscScalar D; /* Damping constant */ 23 PetscScalar Pmax; /* Maximum power output of generator */ 24 PetscScalar PM_min; /* Mean mechanical power input */ 25 PetscScalar lambda; /* correlation time */ 26 PetscScalar q; /* noise strength */ 27 PetscScalar mux; /* Initial average angle */ 28 PetscScalar sigmax; /* Standard deviation of initial angle */ 29 PetscScalar muy; /* Average speed */ 30 PetscScalar sigmay; /* standard deviation of initial speed */ 31 PetscScalar rho; /* Cross-correlation coefficient */ 32 PetscScalar t0; /* Initial time */ 33 PetscScalar tmax; /* Final time */ 34 PetscScalar xmin; /* left boundary of angle */ 35 PetscScalar xmax; /* right boundary of angle */ 36 PetscScalar ymin; /* bottom boundary of speed */ 37 PetscScalar ymax; /* top boundary of speed */ 38 PetscScalar dx; /* x step size */ 39 PetscScalar dy; /* y step size */ 40 PetscInt bc; /* Boundary conditions */ 41 PetscScalar disper_coe; /* Dispersion coefficient */ 42 DM da; 43 } AppCtx; 44 45 PetscErrorCode Parameter_settings(AppCtx *); 46 PetscErrorCode ini_bou(Vec, AppCtx *); 47 PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *); 48 PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *); 49 PetscErrorCode PostStep(TS); 50 51 int main(int argc, char **argv) 52 { 53 Vec x; /* Solution vector */ 54 TS ts; /* Time-stepping context */ 55 AppCtx user; /* Application context */ 56 Mat J; 57 PetscViewer viewer; 58 59 PetscFunctionBeginUser; 60 PetscCall(PetscInitialize(&argc, &argv, "petscopt_ex6", help)); 61 /* Get physics and time parameters */ 62 PetscCall(Parameter_settings(&user)); 63 /* Create a 2D DA with dof = 1 */ 64 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &user.da)); 65 PetscCall(DMSetFromOptions(user.da)); 66 PetscCall(DMSetUp(user.da)); 67 /* Set x and y coordinates */ 68 PetscCall(DMDASetUniformCoordinates(user.da, user.xmin, user.xmax, user.ymin, user.ymax, 0.0, 1.0)); 69 70 /* Get global vector x from DM */ 71 PetscCall(DMCreateGlobalVector(user.da, &x)); 72 73 PetscCall(ini_bou(x, &user)); 74 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "ini_x", FILE_MODE_WRITE, &viewer)); 75 PetscCall(VecView(x, viewer)); 76 PetscCall(PetscViewerDestroy(&viewer)); 77 78 /* Get Jacobian matrix structure from the da */ 79 PetscCall(DMSetMatType(user.da, MATAIJ)); 80 PetscCall(DMCreateMatrix(user.da, &J)); 81 82 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 83 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 84 PetscCall(TSSetIFunction(ts, NULL, IFunction, &user)); 85 PetscCall(TSSetIJacobian(ts, J, J, IJacobian, &user)); 86 PetscCall(TSSetApplicationContext(ts, &user)); 87 PetscCall(TSSetMaxTime(ts, user.tmax)); 88 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 89 PetscCall(TSSetTime(ts, user.t0)); 90 PetscCall(TSSetTimeStep(ts, .005)); 91 PetscCall(TSSetFromOptions(ts)); 92 PetscCall(TSSetPostStep(ts, PostStep)); 93 PetscCall(TSSolve(ts, x)); 94 95 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "fin_x", FILE_MODE_WRITE, &viewer)); 96 PetscCall(VecView(x, viewer)); 97 PetscCall(PetscViewerDestroy(&viewer)); 98 99 PetscCall(VecDestroy(&x)); 100 PetscCall(MatDestroy(&J)); 101 PetscCall(DMDestroy(&user.da)); 102 PetscCall(TSDestroy(&ts)); 103 PetscCall(PetscFinalize()); 104 return 0; 105 } 106 107 PetscErrorCode PostStep(TS ts) 108 { 109 Vec X; 110 AppCtx *user; 111 PetscScalar sum; 112 PetscReal t; 113 114 PetscFunctionBegin; 115 PetscCall(TSGetApplicationContext(ts, &user)); 116 PetscCall(TSGetTime(ts, &t)); 117 PetscCall(TSGetSolution(ts, &X)); 118 PetscCall(VecSum(X, &sum)); 119 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sum(p)*dw*dtheta at t = %3.2f = %3.6f\n", (double)t, (double)(sum * user->dx * user->dy))); 120 PetscFunctionReturn(PETSC_SUCCESS); 121 } 122 123 PetscErrorCode ini_bou(Vec X, AppCtx *user) 124 { 125 DM cda; 126 DMDACoor2d **coors; 127 PetscScalar **p; 128 Vec gc; 129 PetscInt i, j; 130 PetscInt xs, ys, xm, ym, M, N; 131 PetscScalar xi, yi; 132 PetscScalar sigmax = user->sigmax, sigmay = user->sigmay; 133 PetscScalar rho = user->rho; 134 PetscScalar mux = user->mux, muy = user->muy; 135 PetscMPIInt rank; 136 137 PetscFunctionBeginUser; 138 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 139 PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 140 user->dx = (user->xmax - user->xmin) / (M - 1); 141 user->dy = (user->ymax - user->ymin) / (N - 1); 142 PetscCall(DMGetCoordinateDM(user->da, &cda)); 143 PetscCall(DMGetCoordinates(user->da, &gc)); 144 PetscCall(DMDAVecGetArray(cda, gc, &coors)); 145 PetscCall(DMDAVecGetArray(user->da, X, &p)); 146 PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0)); 147 for (i = xs; i < xs + xm; i++) { 148 for (j = ys; j < ys + ym; j++) { 149 xi = coors[j][i].x; 150 yi = coors[j][i].y; 151 if (i == 0 || j == 0 || i == M - 1 || j == N - 1) p[j][i] = 0.0; 152 else 153 p[j][i] = (0.5 / (PETSC_PI * sigmax * sigmay * PetscSqrtScalar(1.0 - rho * rho))) * PetscExpScalar(-0.5 / (1 - rho * rho) * (PetscPowScalar((xi - mux) / sigmax, 2) + PetscPowScalar((yi - muy) / sigmay, 2) - 2 * rho * (xi - mux) * (yi - muy) / (sigmax * sigmay))); 154 } 155 } 156 /* p[N/2+N%2][M/2+M%2] = 1/(user->dx*user->dy); */ 157 158 PetscCall(DMDAVecRestoreArray(cda, gc, &coors)); 159 PetscCall(DMDAVecRestoreArray(user->da, X, &p)); 160 PetscFunctionReturn(PETSC_SUCCESS); 161 } 162 163 /* First advection term */ 164 PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar *p1, AppCtx *user) 165 { 166 PetscScalar f; 167 /* PetscScalar v1,v2,v3,v4,v5,s1,s2,s3; */ 168 PetscFunctionBegin; 169 /* if (i > 2 && i < M-2) { 170 v1 = (y-user->ws)*(p[j][i-2] - p[j][i-3])/user->dx; 171 v2 = (y-user->ws)*(p[j][i-1] - p[j][i-2])/user->dx; 172 v3 = (y-user->ws)*(p[j][i] - p[j][i-1])/user->dx; 173 v4 = (y-user->ws)*(p[j][i+1] - p[j][i])/user->dx; 174 v5 = (y-user->ws)*(p[j][i+1] - p[j][i+2])/user->dx; 175 176 s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3; 177 s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4; 178 s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5; 179 180 *p1 = 0.1*s1 + 0.6*s2 + 0.3*s3; 181 } else *p1 = 0.0; */ 182 f = (y - user->ws); 183 *p1 = f * (p[j][i + 1] - p[j][i - 1]) / (2 * user->dx); 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 /* Second advection term */ 188 PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscInt i, PetscInt j, PetscInt N, PetscScalar *p2, AppCtx *user) 189 { 190 PetscScalar f; 191 /* PetscScalar v1,v2,v3,v4,v5,s1,s2,s3; */ 192 PetscFunctionBegin; 193 /* if (j > 2 && j < N-2) { 194 v1 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-2][i] - p[j-3][i])/user->dy; 195 v2 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-1][i] - p[j-2][i])/user->dy; 196 v3 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j][i] - p[j-1][i])/user->dy; 197 v4 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+1][i] - p[j][i])/user->dy; 198 v5 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+2][i] - p[j+1][i])/user->dy; 199 200 s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3; 201 s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4; 202 s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5; 203 204 *p2 = 0.1*s1 + 0.6*s2 + 0.3*s3; 205 } else *p2 = 0.0; */ 206 f = (user->ws / (2 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(x)); 207 *p2 = f * (p[j + 1][i] - p[j - 1][i]) / (2 * user->dy); 208 PetscFunctionReturn(PETSC_SUCCESS); 209 } 210 211 /* Diffusion term */ 212 PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user) 213 { 214 PetscFunctionBeginUser; 215 216 *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy)); 217 PetscFunctionReturn(PETSC_SUCCESS); 218 } 219 220 PetscErrorCode BoundaryConditions(PetscScalar **p, DMDACoor2d **coors, PetscInt i, PetscInt j, PetscInt M, PetscInt N, PetscScalar **f, AppCtx *user) 221 { 222 PetscScalar fwc, fthetac; 223 PetscScalar w = coors[j][i].y, theta = coors[j][i].x; 224 225 PetscFunctionBeginUser; 226 if (user->bc == 0) { /* Natural boundary condition */ 227 f[j][i] = p[j][i]; 228 } else { /* Steady state boundary condition */ 229 fthetac = user->ws / (2 * user->H) * (user->PM_min - user->Pmax * PetscSinScalar(theta)); 230 fwc = (w * w / 2.0 - user->ws * w); 231 if (i == 0 && j == 0) { /* left bottom corner */ 232 f[j][i] = fwc * (p[j][i + 1] - p[j][i]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j + 1][i] - p[j][i]) / user->dy; 233 } else if (i == 0 && j == N - 1) { /* right bottom corner */ 234 f[j][i] = fwc * (p[j][i + 1] - p[j][i]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j][i] - p[j - 1][i]) / user->dy; 235 } else if (i == M - 1 && j == 0) { /* left top corner */ 236 f[j][i] = fwc * (p[j][i] - p[j][i - 1]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j + 1][i] - p[j][i]) / user->dy; 237 } else if (i == M - 1 && j == N - 1) { /* right top corner */ 238 f[j][i] = fwc * (p[j][i] - p[j][i - 1]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j][i] - p[j - 1][i]) / user->dy; 239 } else if (i == 0) { /* Bottom edge */ 240 f[j][i] = fwc * (p[j][i + 1] - p[j][i]) / (user->dx) + fthetac * p[j][i] - user->disper_coe * (p[j + 1][i] - p[j - 1][i]) / (2 * user->dy); 241 } else if (i == M - 1) { /* Top edge */ 242 f[j][i] = fwc * (p[j][i] - p[j][i - 1]) / (user->dx) + fthetac * p[j][i] - user->disper_coe * (p[j + 1][i] - p[j - 1][i]) / (2 * user->dy); 243 } else if (j == 0) { /* Left edge */ 244 f[j][i] = fwc * (p[j][i + 1] - p[j][i - 1]) / (2 * user->dx) + fthetac * p[j][i] - user->disper_coe * (p[j + 1][i] - p[j][i]) / (user->dy); 245 } else if (j == N - 1) { /* Right edge */ 246 f[j][i] = fwc * (p[j][i + 1] - p[j][i - 1]) / (2 * user->dx) + fthetac * p[j][i] - user->disper_coe * (p[j][i] - p[j - 1][i]) / (user->dy); 247 } 248 } 249 PetscFunctionReturn(PETSC_SUCCESS); 250 } 251 252 PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx) 253 { 254 AppCtx *user = (AppCtx *)ctx; 255 DM cda; 256 DMDACoor2d **coors; 257 PetscScalar **p, **f, **pdot; 258 PetscInt i, j; 259 PetscInt xs, ys, xm, ym, M, N; 260 Vec localX, gc, localXdot; 261 PetscScalar p_adv1, p_adv2, p_diff; 262 263 PetscFunctionBeginUser; 264 PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 265 PetscCall(DMGetCoordinateDM(user->da, &cda)); 266 PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0)); 267 268 PetscCall(DMGetLocalVector(user->da, &localX)); 269 PetscCall(DMGetLocalVector(user->da, &localXdot)); 270 271 PetscCall(DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX)); 272 PetscCall(DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX)); 273 PetscCall(DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot)); 274 PetscCall(DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot)); 275 276 PetscCall(DMGetCoordinatesLocal(user->da, &gc)); 277 278 PetscCall(DMDAVecGetArrayRead(cda, gc, &coors)); 279 PetscCall(DMDAVecGetArrayRead(user->da, localX, &p)); 280 PetscCall(DMDAVecGetArrayRead(user->da, localXdot, &pdot)); 281 PetscCall(DMDAVecGetArray(user->da, F, &f)); 282 283 user->disper_coe = PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda)); 284 for (i = xs; i < xs + xm; i++) { 285 for (j = ys; j < ys + ym; j++) { 286 if (i == 0 || j == 0 || i == M - 1 || j == N - 1) { 287 PetscCall(BoundaryConditions(p, coors, i, j, M, N, f, user)); 288 } else { 289 PetscCall(adv1(p, coors[j][i].y, i, j, M, &p_adv1, user)); 290 PetscCall(adv2(p, coors[j][i].x, i, j, N, &p_adv2, user)); 291 PetscCall(diffuse(p, i, j, t, &p_diff, user)); 292 f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i]; 293 } 294 } 295 } 296 PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &p)); 297 PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &pdot)); 298 PetscCall(DMRestoreLocalVector(user->da, &localX)); 299 PetscCall(DMRestoreLocalVector(user->da, &localXdot)); 300 PetscCall(DMDAVecRestoreArray(user->da, F, &f)); 301 PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors)); 302 303 PetscFunctionReturn(PETSC_SUCCESS); 304 } 305 306 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ctx) 307 { 308 AppCtx *user = (AppCtx *)ctx; 309 DM cda; 310 DMDACoor2d **coors; 311 PetscInt i, j; 312 PetscInt xs, ys, xm, ym, M, N; 313 Vec gc; 314 PetscScalar val[5], xi, yi; 315 MatStencil row, col[5]; 316 PetscScalar c1, c3, c5; 317 318 PetscFunctionBeginUser; 319 PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 320 PetscCall(DMGetCoordinateDM(user->da, &cda)); 321 PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0)); 322 323 PetscCall(DMGetCoordinatesLocal(user->da, &gc)); 324 PetscCall(DMDAVecGetArrayRead(cda, gc, &coors)); 325 for (i = xs; i < xs + xm; i++) { 326 for (j = ys; j < ys + ym; j++) { 327 PetscInt nc = 0; 328 xi = coors[j][i].x; 329 yi = coors[j][i].y; 330 row.i = i; 331 row.j = j; 332 if (i == 0 || j == 0 || i == M - 1 || j == N - 1) { 333 if (user->bc == 0) { 334 col[nc].i = i; 335 col[nc].j = j; 336 val[nc++] = 1.0; 337 } else { 338 PetscScalar fthetac, fwc; 339 fthetac = user->ws / (2 * user->H) * (user->PM_min - user->Pmax * PetscSinScalar(xi)); 340 fwc = (yi * yi / 2.0 - user->ws * yi); 341 if (i == 0 && j == 0) { 342 col[nc].i = i + 1; 343 col[nc].j = j; 344 val[nc++] = fwc / user->dx; 345 col[nc].i = i; 346 col[nc].j = j + 1; 347 val[nc++] = -user->disper_coe / user->dy; 348 col[nc].i = i; 349 col[nc].j = j; 350 val[nc++] = -fwc / user->dx + fthetac + user->disper_coe / user->dy; 351 } else if (i == 0 && j == N - 1) { 352 col[nc].i = i + 1; 353 col[nc].j = j; 354 val[nc++] = fwc / user->dx; 355 col[nc].i = i; 356 col[nc].j = j - 1; 357 val[nc++] = user->disper_coe / user->dy; 358 col[nc].i = i; 359 col[nc].j = j; 360 val[nc++] = -fwc / user->dx + fthetac - user->disper_coe / user->dy; 361 } else if (i == M - 1 && j == 0) { 362 col[nc].i = i - 1; 363 col[nc].j = j; 364 val[nc++] = -fwc / user->dx; 365 col[nc].i = i; 366 col[nc].j = j + 1; 367 val[nc++] = -user->disper_coe / user->dy; 368 col[nc].i = i; 369 col[nc].j = j; 370 val[nc++] = fwc / user->dx + fthetac + user->disper_coe / user->dy; 371 } else if (i == M - 1 && j == N - 1) { 372 col[nc].i = i - 1; 373 col[nc].j = j; 374 val[nc++] = -fwc / user->dx; 375 col[nc].i = i; 376 col[nc].j = j - 1; 377 val[nc++] = user->disper_coe / user->dy; 378 col[nc].i = i; 379 col[nc].j = j; 380 val[nc++] = fwc / user->dx + fthetac - user->disper_coe / user->dy; 381 } else if (i == 0) { 382 col[nc].i = i + 1; 383 col[nc].j = j; 384 val[nc++] = fwc / user->dx; 385 col[nc].i = i; 386 col[nc].j = j + 1; 387 val[nc++] = -user->disper_coe / (2 * user->dy); 388 col[nc].i = i; 389 col[nc].j = j - 1; 390 val[nc++] = user->disper_coe / (2 * user->dy); 391 col[nc].i = i; 392 col[nc].j = j; 393 val[nc++] = -fwc / user->dx + fthetac; 394 } else if (i == M - 1) { 395 col[nc].i = i - 1; 396 col[nc].j = j; 397 val[nc++] = -fwc / user->dx; 398 col[nc].i = i; 399 col[nc].j = j + 1; 400 val[nc++] = -user->disper_coe / (2 * user->dy); 401 col[nc].i = i; 402 col[nc].j = j - 1; 403 val[nc++] = user->disper_coe / (2 * user->dy); 404 col[nc].i = i; 405 col[nc].j = j; 406 val[nc++] = fwc / user->dx + fthetac; 407 } else if (j == 0) { 408 col[nc].i = i + 1; 409 col[nc].j = j; 410 val[nc++] = fwc / (2 * user->dx); 411 col[nc].i = i - 1; 412 col[nc].j = j; 413 val[nc++] = -fwc / (2 * user->dx); 414 col[nc].i = i; 415 col[nc].j = j + 1; 416 val[nc++] = -user->disper_coe / user->dy; 417 col[nc].i = i; 418 col[nc].j = j; 419 val[nc++] = user->disper_coe / user->dy + fthetac; 420 } else if (j == N - 1) { 421 col[nc].i = i + 1; 422 col[nc].j = j; 423 val[nc++] = fwc / (2 * user->dx); 424 col[nc].i = i - 1; 425 col[nc].j = j; 426 val[nc++] = -fwc / (2 * user->dx); 427 col[nc].i = i; 428 col[nc].j = j - 1; 429 val[nc++] = user->disper_coe / user->dy; 430 col[nc].i = i; 431 col[nc].j = j; 432 val[nc++] = -user->disper_coe / user->dy + fthetac; 433 } 434 } 435 } else { 436 c1 = (yi - user->ws) / (2 * user->dx); 437 c3 = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi)) / (2 * user->dy); 438 c5 = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy); 439 col[nc].i = i - 1; 440 col[nc].j = j; 441 val[nc++] = c1; 442 col[nc].i = i + 1; 443 col[nc].j = j; 444 val[nc++] = -c1; 445 col[nc].i = i; 446 col[nc].j = j - 1; 447 val[nc++] = c3 + c5; 448 col[nc].i = i; 449 col[nc].j = j + 1; 450 val[nc++] = -c3 + c5; 451 col[nc].i = i; 452 col[nc].j = j; 453 val[nc++] = -2 * c5 - a; 454 } 455 PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES)); 456 } 457 } 458 PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors)); 459 460 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY)); 461 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY)); 462 if (J != Jpre) { 463 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 464 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 465 } 466 PetscFunctionReturn(PETSC_SUCCESS); 467 } 468 469 PetscErrorCode Parameter_settings(AppCtx *user) 470 { 471 PetscBool flg; 472 473 PetscFunctionBeginUser; 474 475 /* Set default parameters */ 476 user->ws = 1.0; 477 user->H = 5.0; 478 user->Pmax = 2.1; 479 user->PM_min = 1.0; 480 user->lambda = 0.1; 481 user->q = 1.0; 482 user->mux = PetscAsinScalar(user->PM_min / user->Pmax); 483 user->sigmax = 0.1; 484 user->sigmay = 0.1; 485 user->rho = 0.0; 486 user->t0 = 0.0; 487 user->tmax = 2.0; 488 user->xmin = -1.0; 489 user->xmax = 10.0; 490 user->ymin = -1.0; 491 user->ymax = 10.0; 492 user->bc = 0; 493 494 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg)); 495 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg)); 496 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg)); 497 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg)); 498 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg)); 499 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg)); 500 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg)); 501 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmax", &user->sigmax, &flg)); 502 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg)); 503 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmay", &user->sigmay, &flg)); 504 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-rho", &user->rho, &flg)); 505 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-t0", &user->t0, &flg)); 506 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-tmax", &user->tmax, &flg)); 507 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg)); 508 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg)); 509 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg)); 510 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg)); 511 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bc_type", &user->bc, &flg)); 512 user->muy = user->ws; 513 PetscFunctionReturn(PETSC_SUCCESS); 514 } 515 516 /*TEST 517 518 build: 519 requires: !complex 520 521 test: 522 args: -nox -ts_max_steps 2 523 localrunfiles: petscopt_ex6 524 timeoutfactor: 4 525 526 TEST*/ 527