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;

 60:   PetscInitialize(&argc, &argv, "petscopt_ex6", help);
 61:   /* Get physics and time parameters */
 62:   Parameter_settings(&user);
 63:   /* Create a 2D DA with dof = 1 */
 64:   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:   DMSetFromOptions(user.da);
 66:   DMSetUp(user.da);
 67:   /* Set x and y coordinates */
 68:   DMDASetUniformCoordinates(user.da, user.xmin, user.xmax, user.ymin, user.ymax, 0.0, 1.0);

 70:   /* Get global vector x from DM  */
 71:   DMCreateGlobalVector(user.da, &x);

 73:   ini_bou(x, &user);
 74:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "ini_x", FILE_MODE_WRITE, &viewer);
 75:   VecView(x, viewer);
 76:   PetscViewerDestroy(&viewer);

 78:   /* Get Jacobian matrix structure from the da */
 79:   DMSetMatType(user.da, MATAIJ);
 80:   DMCreateMatrix(user.da, &J);

 82:   TSCreate(PETSC_COMM_WORLD, &ts);
 83:   TSSetProblemType(ts, TS_NONLINEAR);
 84:   TSSetIFunction(ts, NULL, IFunction, &user);
 85:   TSSetIJacobian(ts, J, J, IJacobian, &user);
 86:   TSSetApplicationContext(ts, &user);
 87:   TSSetMaxTime(ts, user.tmax);
 88:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
 89:   TSSetTime(ts, user.t0);
 90:   TSSetTimeStep(ts, .005);
 91:   TSSetFromOptions(ts);
 92:   TSSetPostStep(ts, PostStep);
 93:   TSSolve(ts, x);

 95:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "fin_x", FILE_MODE_WRITE, &viewer);
 96:   VecView(x, viewer);
 97:   PetscViewerDestroy(&viewer);

 99:   VecDestroy(&x);
100:   MatDestroy(&J);
101:   DMDestroy(&user.da);
102:   TSDestroy(&ts);
103:   PetscFinalize();
104:   return 0;
105: }

107: PetscErrorCode PostStep(TS ts)
108: {
109:   Vec         X;
110:   AppCtx     *user;
111:   PetscScalar sum;
112:   PetscReal   t;

114:   TSGetApplicationContext(ts, &user);
115:   TSGetTime(ts, &t);
116:   TSGetSolution(ts, &X);
117:   VecSum(X, &sum);
118:   PetscPrintf(PETSC_COMM_WORLD, "sum(p)*dw*dtheta at t = %3.2f = %3.6f\n", (double)t, (double)(sum * user->dx * user->dy));
119:   return 0;
120: }

122: PetscErrorCode ini_bou(Vec X, AppCtx *user)
123: {
124:   DM            cda;
125:   DMDACoor2d  **coors;
126:   PetscScalar **p;
127:   Vec           gc;
128:   PetscInt      i, j;
129:   PetscInt      xs, ys, xm, ym, M, N;
130:   PetscScalar   xi, yi;
131:   PetscScalar   sigmax = user->sigmax, sigmay = user->sigmay;
132:   PetscScalar   rho = user->rho;
133:   PetscScalar   mux = user->mux, muy = user->muy;
134:   PetscMPIInt   rank;

137:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
138:   DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
139:   user->dx = (user->xmax - user->xmin) / (M - 1);
140:   user->dy = (user->ymax - user->ymin) / (N - 1);
141:   DMGetCoordinateDM(user->da, &cda);
142:   DMGetCoordinates(user->da, &gc);
143:   DMDAVecGetArray(cda, gc, &coors);
144:   DMDAVecGetArray(user->da, X, &p);
145:   DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0);
146:   for (i = xs; i < xs + xm; i++) {
147:     for (j = ys; j < ys + ym; j++) {
148:       xi = coors[j][i].x;
149:       yi = coors[j][i].y;
150:       if (i == 0 || j == 0 || i == M - 1 || j == N - 1) p[j][i] = 0.0;
151:       else
152:         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)));
153:     }
154:   }
155:   /*  p[N/2+N%2][M/2+M%2] = 1/(user->dx*user->dy); */

157:   DMDAVecRestoreArray(cda, gc, &coors);
158:   DMDAVecRestoreArray(user->da, X, &p);
159:   return 0;
160: }

162: /* First advection term */
163: PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar *p1, AppCtx *user)
164: {
165:   PetscScalar f;
166:   /*  PetscScalar v1,v2,v3,v4,v5,s1,s2,s3; */
167:   /*  if (i > 2 && i < M-2) {
168:     v1 = (y-user->ws)*(p[j][i-2] - p[j][i-3])/user->dx;
169:     v2 = (y-user->ws)*(p[j][i-1] - p[j][i-2])/user->dx;
170:     v3 = (y-user->ws)*(p[j][i] - p[j][i-1])/user->dx;
171:     v4 = (y-user->ws)*(p[j][i+1] - p[j][i])/user->dx;
172:     v5 = (y-user->ws)*(p[j][i+1] - p[j][i+2])/user->dx;

174:     s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3;
175:     s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4;
176:     s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5;

178:     *p1 = 0.1*s1 + 0.6*s2 + 0.3*s3;
179:     } else *p1 = 0.0; */
180:   f   = (y - user->ws);
181:   *p1 = f * (p[j][i + 1] - p[j][i - 1]) / (2 * user->dx);
182:   return 0;
183: }

185: /* Second advection term */
186: PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscInt i, PetscInt j, PetscInt N, PetscScalar *p2, AppCtx *user)
187: {
188:   PetscScalar f;
189:   /*  PetscScalar v1,v2,v3,v4,v5,s1,s2,s3; */
190:   /*  if (j > 2 && j < N-2) {
191:     v1 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-2][i] - p[j-3][i])/user->dy;
192:     v2 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-1][i] - p[j-2][i])/user->dy;
193:     v3 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j][i] - p[j-1][i])/user->dy;
194:     v4 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+1][i] - p[j][i])/user->dy;
195:     v5 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+2][i] - p[j+1][i])/user->dy;

197:     s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3;
198:     s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4;
199:     s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5;

201:     *p2 = 0.1*s1 + 0.6*s2 + 0.3*s3;
202:     } else *p2 = 0.0; */
203:   f   = (user->ws / (2 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(x));
204:   *p2 = f * (p[j + 1][i] - p[j - 1][i]) / (2 * user->dy);
205:   return 0;
206: }

208: /* Diffusion term */
209: PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user)
210: {

213:   *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy));
214:   return 0;
215: }

217: PetscErrorCode BoundaryConditions(PetscScalar **p, DMDACoor2d **coors, PetscInt i, PetscInt j, PetscInt M, PetscInt N, PetscScalar **f, AppCtx *user)
218: {
219:   PetscScalar fwc, fthetac;
220:   PetscScalar w = coors[j][i].y, theta = coors[j][i].x;

223:   if (user->bc == 0) { /* Natural boundary condition */
224:     f[j][i] = p[j][i];
225:   } else { /* Steady state boundary condition */
226:     fthetac = user->ws / (2 * user->H) * (user->PM_min - user->Pmax * PetscSinScalar(theta));
227:     fwc     = (w * w / 2.0 - user->ws * w);
228:     if (i == 0 && j == 0) { /* left bottom corner */
229:       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;
230:     } else if (i == 0 && j == N - 1) { /* right 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][i] - p[j - 1][i]) / user->dy;
232:     } else if (i == M - 1 && j == 0) { /* left top corner */
233:       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;
234:     } else if (i == M - 1 && j == N - 1) { /* right 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][i] - p[j - 1][i]) / user->dy;
236:     } else if (i == 0) { /* Bottom edge */
237:       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);
238:     } else if (i == M - 1) { /* Top edge */
239:       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);
240:     } else if (j == 0) { /* Left edge */
241:       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);
242:     } else if (j == N - 1) { /* Right 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][i] - p[j - 1][i]) / (user->dy);
244:     }
245:   }
246:   return 0;
247: }

249: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
250: {
251:   AppCtx       *user = (AppCtx *)ctx;
252:   DM            cda;
253:   DMDACoor2d  **coors;
254:   PetscScalar **p, **f, **pdot;
255:   PetscInt      i, j;
256:   PetscInt      xs, ys, xm, ym, M, N;
257:   Vec           localX, gc, localXdot;
258:   PetscScalar   p_adv1, p_adv2, p_diff;

261:   DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
262:   DMGetCoordinateDM(user->da, &cda);
263:   DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0);

265:   DMGetLocalVector(user->da, &localX);
266:   DMGetLocalVector(user->da, &localXdot);

268:   DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX);
269:   DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX);
270:   DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot);
271:   DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot);

273:   DMGetCoordinatesLocal(user->da, &gc);

275:   DMDAVecGetArrayRead(cda, gc, &coors);
276:   DMDAVecGetArrayRead(user->da, localX, &p);
277:   DMDAVecGetArrayRead(user->da, localXdot, &pdot);
278:   DMDAVecGetArray(user->da, F, &f);

280:   user->disper_coe = PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda));
281:   for (i = xs; i < xs + xm; i++) {
282:     for (j = ys; j < ys + ym; j++) {
283:       if (i == 0 || j == 0 || i == M - 1 || j == N - 1) {
284:         BoundaryConditions(p, coors, i, j, M, N, f, user);
285:       } else {
286:         adv1(p, coors[j][i].y, i, j, M, &p_adv1, user);
287:         adv2(p, coors[j][i].x, i, j, N, &p_adv2, user);
288:         diffuse(p, i, j, t, &p_diff, user);
289:         f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
290:       }
291:     }
292:   }
293:   DMDAVecRestoreArrayRead(user->da, localX, &p);
294:   DMDAVecRestoreArrayRead(user->da, localX, &pdot);
295:   DMRestoreLocalVector(user->da, &localX);
296:   DMRestoreLocalVector(user->da, &localXdot);
297:   DMDAVecRestoreArray(user->da, F, &f);
298:   DMDAVecRestoreArrayRead(cda, gc, &coors);

300:   return 0;
301: }

303: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ctx)
304: {
305:   AppCtx      *user = (AppCtx *)ctx;
306:   DM           cda;
307:   DMDACoor2d **coors;
308:   PetscInt     i, j;
309:   PetscInt     xs, ys, xm, ym, M, N;
310:   Vec          gc;
311:   PetscScalar  val[5], xi, yi;
312:   MatStencil   row, col[5];
313:   PetscScalar  c1, c3, c5;

316:   DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
317:   DMGetCoordinateDM(user->da, &cda);
318:   DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0);

320:   DMGetCoordinatesLocal(user->da, &gc);
321:   DMDAVecGetArrayRead(cda, gc, &coors);
322:   for (i = xs; i < xs + xm; i++) {
323:     for (j = ys; j < ys + ym; j++) {
324:       PetscInt nc = 0;
325:       xi          = coors[j][i].x;
326:       yi          = coors[j][i].y;
327:       row.i       = i;
328:       row.j       = j;
329:       if (i == 0 || j == 0 || i == M - 1 || j == N - 1) {
330:         if (user->bc == 0) {
331:           col[nc].i = i;
332:           col[nc].j = j;
333:           val[nc++] = 1.0;
334:         } else {
335:           PetscScalar fthetac, fwc;
336:           fthetac = user->ws / (2 * user->H) * (user->PM_min - user->Pmax * PetscSinScalar(xi));
337:           fwc     = (yi * yi / 2.0 - user->ws * yi);
338:           if (i == 0 && j == 0) {
339:             col[nc].i = i + 1;
340:             col[nc].j = j;
341:             val[nc++] = fwc / user->dx;
342:             col[nc].i = i;
343:             col[nc].j = j + 1;
344:             val[nc++] = -user->disper_coe / user->dy;
345:             col[nc].i = i;
346:             col[nc].j = j;
347:             val[nc++] = -fwc / user->dx + fthetac + user->disper_coe / user->dy;
348:           } else if (i == 0 && j == N - 1) {
349:             col[nc].i = i + 1;
350:             col[nc].j = j;
351:             val[nc++] = fwc / user->dx;
352:             col[nc].i = i;
353:             col[nc].j = j - 1;
354:             val[nc++] = user->disper_coe / user->dy;
355:             col[nc].i = i;
356:             col[nc].j = j;
357:             val[nc++] = -fwc / user->dx + fthetac - user->disper_coe / user->dy;
358:           } else if (i == M - 1 && j == 0) {
359:             col[nc].i = i - 1;
360:             col[nc].j = j;
361:             val[nc++] = -fwc / user->dx;
362:             col[nc].i = i;
363:             col[nc].j = j + 1;
364:             val[nc++] = -user->disper_coe / user->dy;
365:             col[nc].i = i;
366:             col[nc].j = j;
367:             val[nc++] = fwc / user->dx + fthetac + user->disper_coe / user->dy;
368:           } else if (i == M - 1 && j == N - 1) {
369:             col[nc].i = i - 1;
370:             col[nc].j = j;
371:             val[nc++] = -fwc / user->dx;
372:             col[nc].i = i;
373:             col[nc].j = j - 1;
374:             val[nc++] = user->disper_coe / user->dy;
375:             col[nc].i = i;
376:             col[nc].j = j;
377:             val[nc++] = fwc / user->dx + fthetac - user->disper_coe / user->dy;
378:           } else if (i == 0) {
379:             col[nc].i = i + 1;
380:             col[nc].j = j;
381:             val[nc++] = fwc / user->dx;
382:             col[nc].i = i;
383:             col[nc].j = j + 1;
384:             val[nc++] = -user->disper_coe / (2 * user->dy);
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;
390:             val[nc++] = -fwc / user->dx + fthetac;
391:           } else if (i == M - 1) {
392:             col[nc].i = i - 1;
393:             col[nc].j = j;
394:             val[nc++] = -fwc / user->dx;
395:             col[nc].i = i;
396:             col[nc].j = j + 1;
397:             val[nc++] = -user->disper_coe / (2 * user->dy);
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;
403:             val[nc++] = fwc / user->dx + fthetac;
404:           } else if (j == 0) {
405:             col[nc].i = i + 1;
406:             col[nc].j = j;
407:             val[nc++] = fwc / (2 * user->dx);
408:             col[nc].i = i - 1;
409:             col[nc].j = j;
410:             val[nc++] = -fwc / (2 * user->dx);
411:             col[nc].i = i;
412:             col[nc].j = j + 1;
413:             val[nc++] = -user->disper_coe / user->dy;
414:             col[nc].i = i;
415:             col[nc].j = j;
416:             val[nc++] = user->disper_coe / user->dy + fthetac;
417:           } else if (j == N - 1) {
418:             col[nc].i = i + 1;
419:             col[nc].j = j;
420:             val[nc++] = fwc / (2 * user->dx);
421:             col[nc].i = i - 1;
422:             col[nc].j = j;
423:             val[nc++] = -fwc / (2 * user->dx);
424:             col[nc].i = i;
425:             col[nc].j = j - 1;
426:             val[nc++] = user->disper_coe / user->dy;
427:             col[nc].i = i;
428:             col[nc].j = j;
429:             val[nc++] = -user->disper_coe / user->dy + fthetac;
430:           }
431:         }
432:       } else {
433:         c1        = (yi - user->ws) / (2 * user->dx);
434:         c3        = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi)) / (2 * user->dy);
435:         c5        = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy);
436:         col[nc].i = i - 1;
437:         col[nc].j = j;
438:         val[nc++] = c1;
439:         col[nc].i = i + 1;
440:         col[nc].j = j;
441:         val[nc++] = -c1;
442:         col[nc].i = i;
443:         col[nc].j = j - 1;
444:         val[nc++] = c3 + c5;
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;
450:         val[nc++] = -2 * c5 - a;
451:       }
452:       MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES);
453:     }
454:   }
455:   DMDAVecRestoreArrayRead(cda, gc, &coors);

457:   MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
458:   MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
459:   if (J != Jpre) {
460:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
461:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
462:   }
463:   return 0;
464: }

466: PetscErrorCode Parameter_settings(AppCtx *user)
467: {
468:   PetscBool flg;


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;

491:   PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg);
492:   PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg);
493:   PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg);
494:   PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg);
495:   PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg);
496:   PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg);
497:   PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg);
498:   PetscOptionsGetScalar(NULL, NULL, "-sigmax", &user->sigmax, &flg);
499:   PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg);
500:   PetscOptionsGetScalar(NULL, NULL, "-sigmay", &user->sigmay, &flg);
501:   PetscOptionsGetScalar(NULL, NULL, "-rho", &user->rho, &flg);
502:   PetscOptionsGetScalar(NULL, NULL, "-t0", &user->t0, &flg);
503:   PetscOptionsGetScalar(NULL, NULL, "-tmax", &user->tmax, &flg);
504:   PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg);
505:   PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg);
506:   PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg);
507:   PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg);
508:   PetscOptionsGetInt(NULL, NULL, "-bc_type", &user->bc, &flg);
509:   user->muy = user->ws;
510:   return 0;
511: }

513: /*TEST

515:    build:
516:       requires: !complex

518:    test:
519:       args: -nox -ts_max_steps 2
520:       localrunfiles: petscopt_ex6
521:       timeoutfactor: 4

523: TEST*/