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
main(int argc,char ** argv)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 nonzero 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
PostStep(TS ts)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
ini_bou(Vec X,AppCtx * user)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 */
adv1(PetscScalar ** p,PetscScalar y,PetscInt i,PetscInt j,PetscInt M,PetscScalar * p1,AppCtx * user)164 PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar *p1, AppCtx *user)
165 {
166 PetscScalar f; /* v1,v2,v3,v4,v5,s1,s2,s3; */
167
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 */
adv2(PetscScalar ** p,PetscScalar x,PetscInt i,PetscInt j,PetscInt N,PetscScalar * p2,AppCtx * user)188 PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscInt i, PetscInt j, PetscInt N, PetscScalar *p2, AppCtx *user)
189 {
190 PetscScalar f; /* v1,v2,v3,v4,v5,s1,s2,s3; */
191
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 */
diffuse(PetscScalar ** p,PetscInt i,PetscInt j,PetscReal t,PetscScalar * p_diff,AppCtx * user)212 PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user)
213 {
214 PetscFunctionBeginUser;
215 *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy));
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
BoundaryConditions(PetscScalar ** p,DMDACoor2d ** coors,PetscInt i,PetscInt j,PetscInt M,PetscInt N,PetscScalar ** f,AppCtx * user)219 PetscErrorCode BoundaryConditions(PetscScalar **p, DMDACoor2d **coors, PetscInt i, PetscInt j, PetscInt M, PetscInt N, PetscScalar **f, AppCtx *user)
220 {
221 PetscScalar fwc, fthetac;
222 PetscScalar w = coors[j][i].y, theta = coors[j][i].x;
223
224 PetscFunctionBeginUser;
225 if (user->bc == 0) { /* Natural boundary condition */
226 f[j][i] = p[j][i];
227 } else { /* Steady state boundary condition */
228 fthetac = user->ws / (2 * user->H) * (user->PM_min - user->Pmax * PetscSinScalar(theta));
229 fwc = (w * w / 2.0 - user->ws * w);
230 if (i == 0 && j == 0) { /* left bottom corner */
231 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;
232 } else if (i == 0 && j == N - 1) { /* right bottom corner */
233 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;
234 } else if (i == M - 1 && j == 0) { /* left top corner */
235 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;
236 } else if (i == M - 1 && j == N - 1) { /* right top corner */
237 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;
238 } else if (i == 0) { /* Bottom edge */
239 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);
240 } else if (i == M - 1) { /* Top edge */
241 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);
242 } else if (j == 0) { /* Left edge */
243 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);
244 } else if (j == N - 1) { /* Right edge */
245 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);
246 }
247 }
248 PetscFunctionReturn(PETSC_SUCCESS);
249 }
250
IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,PetscCtx ctx)251 PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx)
252 {
253 AppCtx *user = (AppCtx *)ctx;
254 DM cda;
255 DMDACoor2d **coors;
256 PetscScalar **p, **f, **pdot;
257 PetscInt i, j;
258 PetscInt xs, ys, xm, ym, M, N;
259 Vec localX, gc, localXdot;
260 PetscScalar p_adv1, p_adv2, p_diff;
261
262 PetscFunctionBeginUser;
263 PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
264 PetscCall(DMGetCoordinateDM(user->da, &cda));
265 PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
266
267 PetscCall(DMGetLocalVector(user->da, &localX));
268 PetscCall(DMGetLocalVector(user->da, &localXdot));
269
270 PetscCall(DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX));
271 PetscCall(DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX));
272 PetscCall(DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot));
273 PetscCall(DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot));
274
275 PetscCall(DMGetCoordinatesLocal(user->da, &gc));
276
277 PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
278 PetscCall(DMDAVecGetArrayRead(user->da, localX, &p));
279 PetscCall(DMDAVecGetArrayRead(user->da, localXdot, &pdot));
280 PetscCall(DMDAVecGetArray(user->da, F, &f));
281
282 user->disper_coe = PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda));
283 for (i = xs; i < xs + xm; i++) {
284 for (j = ys; j < ys + ym; j++) {
285 if (i == 0 || j == 0 || i == M - 1 || j == N - 1) {
286 PetscCall(BoundaryConditions(p, coors, i, j, M, N, f, user));
287 } else {
288 PetscCall(adv1(p, coors[j][i].y, i, j, M, &p_adv1, user));
289 PetscCall(adv2(p, coors[j][i].x, i, j, N, &p_adv2, user));
290 PetscCall(diffuse(p, i, j, t, &p_diff, user));
291 f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
292 }
293 }
294 }
295 PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &p));
296 PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &pdot));
297 PetscCall(DMRestoreLocalVector(user->da, &localX));
298 PetscCall(DMRestoreLocalVector(user->da, &localXdot));
299 PetscCall(DMDAVecRestoreArray(user->da, F, &f));
300 PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
301 PetscFunctionReturn(PETSC_SUCCESS);
302 }
303
IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,PetscCtx ctx)304 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, PetscCtx ctx)
305 {
306 AppCtx *user = (AppCtx *)ctx;
307 DM cda;
308 DMDACoor2d **coors;
309 PetscInt i, j;
310 PetscInt xs, ys, xm, ym, M, N;
311 Vec gc;
312 PetscScalar val[5], xi, yi;
313 MatStencil row, col[5];
314 PetscScalar c1, c3, c5;
315
316 PetscFunctionBeginUser;
317 PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
318 PetscCall(DMGetCoordinateDM(user->da, &cda));
319 PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
320
321 PetscCall(DMGetCoordinatesLocal(user->da, &gc));
322 PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
323 for (i = xs; i < xs + xm; i++) {
324 for (j = ys; j < ys + ym; j++) {
325 PetscInt nc = 0;
326 xi = coors[j][i].x;
327 yi = coors[j][i].y;
328 row.i = i;
329 row.j = j;
330 if (i == 0 || j == 0 || i == M - 1 || j == N - 1) {
331 if (user->bc == 0) {
332 col[nc].i = i;
333 col[nc].j = j;
334 val[nc++] = 1.0;
335 } else {
336 PetscScalar fthetac, fwc;
337 fthetac = user->ws / (2 * user->H) * (user->PM_min - user->Pmax * PetscSinScalar(xi));
338 fwc = (yi * yi / 2.0 - user->ws * yi);
339 if (i == 0 && j == 0) {
340 col[nc].i = i + 1;
341 col[nc].j = j;
342 val[nc++] = fwc / user->dx;
343 col[nc].i = i;
344 col[nc].j = j + 1;
345 val[nc++] = -user->disper_coe / user->dy;
346 col[nc].i = i;
347 col[nc].j = j;
348 val[nc++] = -fwc / user->dx + fthetac + user->disper_coe / user->dy;
349 } else if (i == 0 && j == N - 1) {
350 col[nc].i = i + 1;
351 col[nc].j = j;
352 val[nc++] = fwc / user->dx;
353 col[nc].i = i;
354 col[nc].j = j - 1;
355 val[nc++] = user->disper_coe / user->dy;
356 col[nc].i = i;
357 col[nc].j = j;
358 val[nc++] = -fwc / user->dx + fthetac - user->disper_coe / user->dy;
359 } else if (i == M - 1 && j == 0) {
360 col[nc].i = i - 1;
361 col[nc].j = j;
362 val[nc++] = -fwc / user->dx;
363 col[nc].i = i;
364 col[nc].j = j + 1;
365 val[nc++] = -user->disper_coe / user->dy;
366 col[nc].i = i;
367 col[nc].j = j;
368 val[nc++] = fwc / user->dx + fthetac + user->disper_coe / user->dy;
369 } else if (i == M - 1 && j == N - 1) {
370 col[nc].i = i - 1;
371 col[nc].j = j;
372 val[nc++] = -fwc / user->dx;
373 col[nc].i = i;
374 col[nc].j = j - 1;
375 val[nc++] = user->disper_coe / user->dy;
376 col[nc].i = i;
377 col[nc].j = j;
378 val[nc++] = fwc / user->dx + fthetac - user->disper_coe / user->dy;
379 } else if (i == 0) {
380 col[nc].i = i + 1;
381 col[nc].j = j;
382 val[nc++] = fwc / user->dx;
383 col[nc].i = i;
384 col[nc].j = j + 1;
385 val[nc++] = -user->disper_coe / (2 * user->dy);
386 col[nc].i = i;
387 col[nc].j = j - 1;
388 val[nc++] = user->disper_coe / (2 * user->dy);
389 col[nc].i = i;
390 col[nc].j = j;
391 val[nc++] = -fwc / user->dx + fthetac;
392 } else if (i == M - 1) {
393 col[nc].i = i - 1;
394 col[nc].j = j;
395 val[nc++] = -fwc / user->dx;
396 col[nc].i = i;
397 col[nc].j = j + 1;
398 val[nc++] = -user->disper_coe / (2 * user->dy);
399 col[nc].i = i;
400 col[nc].j = j - 1;
401 val[nc++] = user->disper_coe / (2 * user->dy);
402 col[nc].i = i;
403 col[nc].j = j;
404 val[nc++] = fwc / user->dx + fthetac;
405 } else if (j == 0) {
406 col[nc].i = i + 1;
407 col[nc].j = j;
408 val[nc++] = fwc / (2 * user->dx);
409 col[nc].i = i - 1;
410 col[nc].j = j;
411 val[nc++] = -fwc / (2 * user->dx);
412 col[nc].i = i;
413 col[nc].j = j + 1;
414 val[nc++] = -user->disper_coe / user->dy;
415 col[nc].i = i;
416 col[nc].j = j;
417 val[nc++] = user->disper_coe / user->dy + fthetac;
418 } else if (j == N - 1) {
419 col[nc].i = i + 1;
420 col[nc].j = j;
421 val[nc++] = fwc / (2 * user->dx);
422 col[nc].i = i - 1;
423 col[nc].j = j;
424 val[nc++] = -fwc / (2 * user->dx);
425 col[nc].i = i;
426 col[nc].j = j - 1;
427 val[nc++] = user->disper_coe / user->dy;
428 col[nc].i = i;
429 col[nc].j = j;
430 val[nc++] = -user->disper_coe / user->dy + fthetac;
431 }
432 }
433 } else {
434 c1 = (yi - user->ws) / (2 * user->dx);
435 c3 = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi)) / (2 * user->dy);
436 c5 = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy);
437 col[nc].i = i - 1;
438 col[nc].j = j;
439 val[nc++] = c1;
440 col[nc].i = i + 1;
441 col[nc].j = j;
442 val[nc++] = -c1;
443 col[nc].i = i;
444 col[nc].j = j - 1;
445 val[nc++] = c3 + c5;
446 col[nc].i = i;
447 col[nc].j = j + 1;
448 val[nc++] = -c3 + c5;
449 col[nc].i = i;
450 col[nc].j = j;
451 val[nc++] = -2 * c5 - a;
452 }
453 PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES));
454 }
455 }
456 PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
457
458 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
459 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
460 if (J != Jpre) {
461 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
462 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
463 }
464 PetscFunctionReturn(PETSC_SUCCESS);
465 }
466
Parameter_settings(AppCtx * user)467 PetscErrorCode Parameter_settings(AppCtx *user)
468 {
469 PetscBool flg;
470
471 PetscFunctionBeginUser;
472 /* Set default parameters */
473 user->ws = 1.0;
474 user->H = 5.0;
475 user->Pmax = 2.1;
476 user->PM_min = 1.0;
477 user->lambda = 0.1;
478 user->q = 1.0;
479 user->mux = PetscAsinScalar(user->PM_min / user->Pmax);
480 user->sigmax = 0.1;
481 user->sigmay = 0.1;
482 user->rho = 0.0;
483 user->t0 = 0.0;
484 user->tmax = 2.0;
485 user->xmin = -1.0;
486 user->xmax = 10.0;
487 user->ymin = -1.0;
488 user->ymax = 10.0;
489 user->bc = 0;
490
491 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg));
492 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg));
493 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg));
494 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg));
495 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg));
496 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg));
497 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg));
498 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmax", &user->sigmax, &flg));
499 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg));
500 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmay", &user->sigmay, &flg));
501 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-rho", &user->rho, &flg));
502 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-t0", &user->t0, &flg));
503 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-tmax", &user->tmax, &flg));
504 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg));
505 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg));
506 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg));
507 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg));
508 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bc_type", &user->bc, &flg));
509 user->muy = user->ws;
510 PetscFunctionReturn(PETSC_SUCCESS);
511 }
512
513 /*TEST
514
515 build:
516 requires: !complex
517
518 test:
519 args: -nox -ts_max_steps 2
520 localrunfiles: petscopt_ex6
521 timeoutfactor: 4
522
523 TEST*/
524