Lines Matching +full:- +full:j
1 static char help[] = "Time-dependent PDE in 2d for calculating joint PDF. \n";
3 p_t = -x_t*p_x -y_t*p_y + f(t)*p_yy
5 x_t = (y - ws) y_t = (ws/2H)*(Pm - Pmax*sin(x))
7 Boundary conditions: -bc_type 0 => Zero dirichlet boundary
8 -bc_type 1 => Steady state boundary condition
17 User-defined data structures and routines
31 PetscScalar rho; /* Cross-correlation coefficient */
54 TS ts; /* Time-stepping context */ in main()
56 Mat J; in main() local
80 PetscCall(DMCreateMatrix(user.da, &J)); in main()
85 PetscCall(TSSetIJacobian(ts, J, J, IJacobian, &user)); in main()
100 PetscCall(MatDestroy(&J)); in main()
119 …WORLD, "sum(p)*dw*dtheta at t = %3.2f = %3.6f\n", (double)t, (double)(sum * user->dx * user->dy))); in PostStep()
129 PetscInt i, j; in ini_bou() local
132 PetscScalar sigmax = user->sigmax, sigmay = user->sigmay; in ini_bou()
133 PetscScalar rho = user->rho; in ini_bou()
134 PetscScalar mux = user->mux, muy = user->muy; in ini_bou()
139 …PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL… in ini_bou()
140 user->dx = (user->xmax - user->xmin) / (M - 1); in ini_bou()
141 user->dy = (user->ymax - user->ymin) / (N - 1); in ini_bou()
142 PetscCall(DMGetCoordinateDM(user->da, &cda)); in ini_bou()
143 PetscCall(DMGetCoordinates(user->da, &gc)); in ini_bou()
145 PetscCall(DMDAVecGetArray(user->da, X, &p)); in ini_bou()
148 for (j = ys; j < ys + ym; j++) { in ini_bou()
149 xi = coors[j][i].x; in ini_bou()
150 yi = coors[j][i].y; in ini_bou()
151 if (i == 0 || j == 0 || i == M - 1 || j == N - 1) p[j][i] = 0.0; in ini_bou()
153 …j][i] = (0.5 / (PETSC_PI * sigmax * sigmay * PetscSqrtScalar(1.0 - rho * rho))) * PetscExpScalar(-… in ini_bou()
156 /* p[N/2+N%2][M/2+M%2] = 1/(user->dx*user->dy); */ in ini_bou()
159 PetscCall(DMDAVecRestoreArray(user->da, X, &p)); in ini_bou()
164 PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar… in adv1() argument
169 /* if (i > 2 && i < M-2) { in adv1()
170 v1 = (y-user->ws)*(p[j][i-2] - p[j][i-3])/user->dx; in adv1()
171 v2 = (y-user->ws)*(p[j][i-1] - p[j][i-2])/user->dx; in adv1()
172 v3 = (y-user->ws)*(p[j][i] - p[j][i-1])/user->dx; in adv1()
173 v4 = (y-user->ws)*(p[j][i+1] - p[j][i])/user->dx; in adv1()
174 v5 = (y-user->ws)*(p[j][i+1] - p[j][i+2])/user->dx; in adv1()
176 s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3; in adv1()
177 s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4; in adv1()
178 s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5; in adv1()
182 f = (y - user->ws); in adv1()
183 *p1 = f * (p[j][i + 1] - p[j][i - 1]) / (2 * user->dx); in adv1()
188 PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscInt i, PetscInt j, PetscInt N, PetscScalar… in adv2() argument
193 /* if (j > 2 && j < N-2) { in adv2()
194 v1 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-2][i] - p[j-3][i])/user->dy; in adv2()
195 v2 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-1][i] - p[j-2][i])/user->dy; in adv2()
196 v3 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j][i] - p[j-1][i])/user->dy; in adv2()
197 v4 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+1][i] - p[j][i])/user->dy; in adv2()
198 v5 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+2][i] - p[j+1][i])/user->dy; in adv2()
200 s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3; in adv2()
201 s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4; in adv2()
202 s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5; in adv2()
206 f = (user->ws / (2 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(x)); in adv2()
207 *p2 = f * (p[j + 1][i] - p[j - 1][i]) / (2 * user->dy); in adv2()
212 PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, A… in diffuse() argument
215 *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy)); in diffuse()
219 PetscErrorCode BoundaryConditions(PetscScalar **p, DMDACoor2d **coors, PetscInt i, PetscInt j, Pets… in BoundaryConditions() argument
222 PetscScalar w = coors[j][i].y, theta = coors[j][i].x; in BoundaryConditions()
225 if (user->bc == 0) { /* Natural boundary condition */ in BoundaryConditions()
226 f[j][i] = p[j][i]; in BoundaryConditions()
228 fthetac = user->ws / (2 * user->H) * (user->PM_min - user->Pmax * PetscSinScalar(theta)); in BoundaryConditions()
229 fwc = (w * w / 2.0 - user->ws * w); in BoundaryConditions()
230 if (i == 0 && j == 0) { /* left bottom corner */ in BoundaryConditions()
231 …f[j][i] = fwc * (p[j][i + 1] - p[j][i]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j +… in BoundaryConditions()
232 } else if (i == 0 && j == N - 1) { /* right bottom corner */ in BoundaryConditions()
233 …f[j][i] = fwc * (p[j][i + 1] - p[j][i]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j][… in BoundaryConditions()
234 } else if (i == M - 1 && j == 0) { /* left top corner */ in BoundaryConditions()
235 …f[j][i] = fwc * (p[j][i] - p[j][i - 1]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j +… in BoundaryConditions()
236 } else if (i == M - 1 && j == N - 1) { /* right top corner */ in BoundaryConditions()
237 …f[j][i] = fwc * (p[j][i] - p[j][i - 1]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j][… in BoundaryConditions()
239 …f[j][i] = fwc * (p[j][i + 1] - p[j][i]) / (user->dx) + fthetac * p[j][i] - user->disper_coe * (p[j… in BoundaryConditions()
240 } else if (i == M - 1) { /* Top edge */ in BoundaryConditions()
241 …f[j][i] = fwc * (p[j][i] - p[j][i - 1]) / (user->dx) + fthetac * p[j][i] - user->disper_coe * (p[j… in BoundaryConditions()
242 } else if (j == 0) { /* Left edge */ in BoundaryConditions()
243 …f[j][i] = fwc * (p[j][i + 1] - p[j][i - 1]) / (2 * user->dx) + fthetac * p[j][i] - user->disper_co… in BoundaryConditions()
244 } else if (j == N - 1) { /* Right edge */ in BoundaryConditions()
245 …f[j][i] = fwc * (p[j][i + 1] - p[j][i - 1]) / (2 * user->dx) + fthetac * p[j][i] - user->disper_co… in BoundaryConditions()
257 PetscInt i, j; in IFunction() local
263 …PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL… in IFunction()
264 PetscCall(DMGetCoordinateDM(user->da, &cda)); in IFunction()
267 PetscCall(DMGetLocalVector(user->da, &localX)); in IFunction()
268 PetscCall(DMGetLocalVector(user->da, &localXdot)); in IFunction()
270 PetscCall(DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX)); in IFunction()
271 PetscCall(DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX)); in IFunction()
272 PetscCall(DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot)); in IFunction()
273 PetscCall(DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot)); in IFunction()
275 PetscCall(DMGetCoordinatesLocal(user->da, &gc)); in IFunction()
278 PetscCall(DMDAVecGetArrayRead(user->da, localX, &p)); in IFunction()
279 PetscCall(DMDAVecGetArrayRead(user->da, localXdot, &pdot)); in IFunction()
280 PetscCall(DMDAVecGetArray(user->da, F, &f)); in IFunction()
282 …user->disper_coe = PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 -… in IFunction()
284 for (j = ys; j < ys + ym; j++) { in IFunction()
285 if (i == 0 || j == 0 || i == M - 1 || j == N - 1) { in IFunction()
286 PetscCall(BoundaryConditions(p, coors, i, j, M, N, f, user)); in IFunction()
288 PetscCall(adv1(p, coors[j][i].y, i, j, M, &p_adv1, user)); in IFunction()
289 PetscCall(adv2(p, coors[j][i].x, i, j, N, &p_adv2, user)); in IFunction()
290 PetscCall(diffuse(p, i, j, t, &p_diff, user)); in IFunction()
291 f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i]; in IFunction()
295 PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &p)); in IFunction()
296 PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &pdot)); in IFunction()
297 PetscCall(DMRestoreLocalVector(user->da, &localX)); in IFunction()
298 PetscCall(DMRestoreLocalVector(user->da, &localXdot)); in IFunction()
299 PetscCall(DMDAVecRestoreArray(user->da, F, &f)); in IFunction()
304 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, PetscCt… in IJacobian() argument
309 PetscInt i, j; in IJacobian() local
317 …PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL… in IJacobian()
318 PetscCall(DMGetCoordinateDM(user->da, &cda)); in IJacobian()
321 PetscCall(DMGetCoordinatesLocal(user->da, &gc)); in IJacobian()
324 for (j = ys; j < ys + ym; j++) { in IJacobian()
326 xi = coors[j][i].x; in IJacobian()
327 yi = coors[j][i].y; in IJacobian()
329 row.j = j; in IJacobian()
330 if (i == 0 || j == 0 || i == M - 1 || j == N - 1) { in IJacobian()
331 if (user->bc == 0) { in IJacobian()
333 col[nc].j = j; in IJacobian()
337 fthetac = user->ws / (2 * user->H) * (user->PM_min - user->Pmax * PetscSinScalar(xi)); in IJacobian()
338 fwc = (yi * yi / 2.0 - user->ws * yi); in IJacobian()
339 if (i == 0 && j == 0) { in IJacobian()
341 col[nc].j = j; in IJacobian()
342 val[nc++] = fwc / user->dx; in IJacobian()
344 col[nc].j = j + 1; in IJacobian()
345 val[nc++] = -user->disper_coe / user->dy; in IJacobian()
347 col[nc].j = j; in IJacobian()
348 val[nc++] = -fwc / user->dx + fthetac + user->disper_coe / user->dy; in IJacobian()
349 } else if (i == 0 && j == N - 1) { in IJacobian()
351 col[nc].j = j; in IJacobian()
352 val[nc++] = fwc / user->dx; in IJacobian()
354 col[nc].j = j - 1; in IJacobian()
355 val[nc++] = user->disper_coe / user->dy; in IJacobian()
357 col[nc].j = j; in IJacobian()
358 val[nc++] = -fwc / user->dx + fthetac - user->disper_coe / user->dy; in IJacobian()
359 } else if (i == M - 1 && j == 0) { in IJacobian()
360 col[nc].i = i - 1; in IJacobian()
361 col[nc].j = j; in IJacobian()
362 val[nc++] = -fwc / user->dx; in IJacobian()
364 col[nc].j = j + 1; in IJacobian()
365 val[nc++] = -user->disper_coe / user->dy; in IJacobian()
367 col[nc].j = j; in IJacobian()
368 val[nc++] = fwc / user->dx + fthetac + user->disper_coe / user->dy; in IJacobian()
369 } else if (i == M - 1 && j == N - 1) { in IJacobian()
370 col[nc].i = i - 1; in IJacobian()
371 col[nc].j = j; in IJacobian()
372 val[nc++] = -fwc / user->dx; in IJacobian()
374 col[nc].j = j - 1; in IJacobian()
375 val[nc++] = user->disper_coe / user->dy; in IJacobian()
377 col[nc].j = j; in IJacobian()
378 val[nc++] = fwc / user->dx + fthetac - user->disper_coe / user->dy; in IJacobian()
381 col[nc].j = j; in IJacobian()
382 val[nc++] = fwc / user->dx; in IJacobian()
384 col[nc].j = j + 1; in IJacobian()
385 val[nc++] = -user->disper_coe / (2 * user->dy); in IJacobian()
387 col[nc].j = j - 1; in IJacobian()
388 val[nc++] = user->disper_coe / (2 * user->dy); in IJacobian()
390 col[nc].j = j; in IJacobian()
391 val[nc++] = -fwc / user->dx + fthetac; in IJacobian()
392 } else if (i == M - 1) { in IJacobian()
393 col[nc].i = i - 1; in IJacobian()
394 col[nc].j = j; in IJacobian()
395 val[nc++] = -fwc / user->dx; in IJacobian()
397 col[nc].j = j + 1; in IJacobian()
398 val[nc++] = -user->disper_coe / (2 * user->dy); in IJacobian()
400 col[nc].j = j - 1; in IJacobian()
401 val[nc++] = user->disper_coe / (2 * user->dy); in IJacobian()
403 col[nc].j = j; in IJacobian()
404 val[nc++] = fwc / user->dx + fthetac; in IJacobian()
405 } else if (j == 0) { in IJacobian()
407 col[nc].j = j; in IJacobian()
408 val[nc++] = fwc / (2 * user->dx); in IJacobian()
409 col[nc].i = i - 1; in IJacobian()
410 col[nc].j = j; in IJacobian()
411 val[nc++] = -fwc / (2 * user->dx); in IJacobian()
413 col[nc].j = j + 1; in IJacobian()
414 val[nc++] = -user->disper_coe / user->dy; in IJacobian()
416 col[nc].j = j; in IJacobian()
417 val[nc++] = user->disper_coe / user->dy + fthetac; in IJacobian()
418 } else if (j == N - 1) { in IJacobian()
420 col[nc].j = j; in IJacobian()
421 val[nc++] = fwc / (2 * user->dx); in IJacobian()
422 col[nc].i = i - 1; in IJacobian()
423 col[nc].j = j; in IJacobian()
424 val[nc++] = -fwc / (2 * user->dx); in IJacobian()
426 col[nc].j = j - 1; in IJacobian()
427 val[nc++] = user->disper_coe / user->dy; in IJacobian()
429 col[nc].j = j; in IJacobian()
430 val[nc++] = -user->disper_coe / user->dy + fthetac; in IJacobian()
434 c1 = (yi - user->ws) / (2 * user->dx); in IJacobian()
435 …c3 = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi)) / (2 *… in IJacobian()
436 …etscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t /… in IJacobian()
437 col[nc].i = i - 1; in IJacobian()
438 col[nc].j = j; in IJacobian()
441 col[nc].j = j; in IJacobian()
442 val[nc++] = -c1; in IJacobian()
444 col[nc].j = j - 1; in IJacobian()
447 col[nc].j = j + 1; in IJacobian()
448 val[nc++] = -c3 + c5; in IJacobian()
450 col[nc].j = j; in IJacobian()
451 val[nc++] = -2 * c5 - a; in IJacobian()
460 if (J != Jpre) { in IJacobian()
461 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); in IJacobian()
462 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); in IJacobian()
473 user->ws = 1.0; in Parameter_settings()
474 user->H = 5.0; in Parameter_settings()
475 user->Pmax = 2.1; in Parameter_settings()
476 user->PM_min = 1.0; in Parameter_settings()
477 user->lambda = 0.1; in Parameter_settings()
478 user->q = 1.0; in Parameter_settings()
479 user->mux = PetscAsinScalar(user->PM_min / user->Pmax); in Parameter_settings()
480 user->sigmax = 0.1; in Parameter_settings()
481 user->sigmay = 0.1; in Parameter_settings()
482 user->rho = 0.0; in Parameter_settings()
483 user->t0 = 0.0; in Parameter_settings()
484 user->tmax = 2.0; in Parameter_settings()
485 user->xmin = -1.0; in Parameter_settings()
486 user->xmax = 10.0; in Parameter_settings()
487 user->ymin = -1.0; in Parameter_settings()
488 user->ymax = 10.0; in Parameter_settings()
489 user->bc = 0; in Parameter_settings()
491 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg)); in Parameter_settings()
492 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg)); in Parameter_settings()
493 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg)); in Parameter_settings()
494 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg)); in Parameter_settings()
495 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg)); in Parameter_settings()
496 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg)); in Parameter_settings()
497 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg)); in Parameter_settings()
498 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmax", &user->sigmax, &flg)); in Parameter_settings()
499 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg)); in Parameter_settings()
500 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmay", &user->sigmay, &flg)); in Parameter_settings()
501 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-rho", &user->rho, &flg)); in Parameter_settings()
502 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-t0", &user->t0, &flg)); in Parameter_settings()
503 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-tmax", &user->tmax, &flg)); in Parameter_settings()
504 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg)); in Parameter_settings()
505 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg)); in Parameter_settings()
506 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg)); in Parameter_settings()
507 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg)); in Parameter_settings()
508 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bc_type", &user->bc, &flg)); in Parameter_settings()
509 user->muy = user->ws; in Parameter_settings()
519 args: -nox -ts_max_steps 2