Actual source code: ex6.c
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))
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: */
12: #include <petscdm.h>
13: #include <petscdmda.h>
14: #include <petscts.h>
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;
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);
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;
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));
70: /* Get global vector x from DM */
71: PetscCall(DMCreateGlobalVector(user.da, &x));
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));
78: /* Get Jacobian matrix structure from the da */
79: PetscCall(DMSetMatType(user.da, MATAIJ));
80: PetscCall(DMCreateMatrix(user.da, &J));
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));
95: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "fin_x", FILE_MODE_WRITE, &viewer));
96: PetscCall(VecView(x, viewer));
97: PetscCall(PetscViewerDestroy(&viewer));
99: PetscCall(VecDestroy(&x));
100: PetscCall(MatDestroy(&J));
101: PetscCall(DMDestroy(&user.da));
102: PetscCall(TSDestroy(&ts));
103: PetscCall(PetscFinalize());
104: return 0;
105: }
107: PetscErrorCode PostStep(TS ts)
108: {
109: Vec X;
110: AppCtx *user;
111: PetscScalar sum;
112: PetscReal t;
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: }
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;
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); */
158: PetscCall(DMDAVecRestoreArray(cda, gc, &coors));
159: PetscCall(DMDAVecRestoreArray(user->da, X, &p));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
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;
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;
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: }
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;
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;
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: }
211: /* Diffusion term */
212: PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user)
213: {
214: PetscFunctionBeginUser;
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: }
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;
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: }
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;
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));
268: PetscCall(DMGetLocalVector(user->da, &localX));
269: PetscCall(DMGetLocalVector(user->da, &localXdot));
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));
276: PetscCall(DMGetCoordinatesLocal(user->da, &gc));
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));
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));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
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;
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));
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));
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: }
469: PetscErrorCode Parameter_settings(AppCtx *user)
470: {
471: PetscBool flg;
473: PetscFunctionBeginUser;
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;
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: }
516: /*TEST
518: build:
519: requires: !complex
521: test:
522: args: -nox -ts_max_steps 2
523: localrunfiles: petscopt_ex6
524: timeoutfactor: 4
526: TEST*/