| 1 |
REQUIRE "ivpsystem.a4l"; |
| 2 |
REQUIRE "atoms.a4l"; |
| 3 |
|
| 4 |
(* |
| 5 |
* ----------------------------------------------------------------- |
| 6 |
* Example problem for IDA: 2D heat equation, serial, GMRES. |
| 7 |
* |
| 8 |
* This example solves a discretized 2D heat equation problem. |
| 9 |
* This version uses the Krylov solver IDASpgmr. |
| 10 |
* |
| 11 |
* The DAE system solved is a spatial discretization of the PDE |
| 12 |
* du/dt = d^2u/dx^2 + d^2u/dy^2 |
| 13 |
* on the unit square. The boundary condition is u = 0 on all edges. |
| 14 |
* Initial conditions are given by u = 16 x (1 - x) y (1 - y). The |
| 15 |
* PDE is treated with central differences on a uniform M x M grid. |
| 16 |
* The values of u at the interior points satisfy ODEs, and |
| 17 |
* equations u = 0 at the boundaries are appended, to form a DAE |
| 18 |
* system of size N = M^2. Here M = 10. |
| 19 |
* |
| 20 |
* The system is solved with IDA/IDAS using the Krylov linear solver |
| 21 |
* IDASPGMR. The preconditioner uses the diagonal elements of the |
| 22 |
* Jacobian only. Routines for preconditioning, required by |
| 23 |
* IDASPGMR, are supplied here. The constraints u >= 0 are posed |
| 24 |
* for all components. Output is taken at t = 0, .01, .02, .04, |
| 25 |
* ..., 10.24. Two cases are run -- with the Gram-Schmidt type |
| 26 |
* being Modified in the first case, and Classical in the second. |
| 27 |
* The second run uses IDAReInit and IDAReInitSpgmr. |
| 28 |
* ----------------------------------------------------------------- |
| 29 |
*) |
| 30 |
MODEL idakryx; |
| 31 |
|
| 32 |
M IS_A integer_constant; |
| 33 |
M :== 5; |
| 34 |
|
| 35 |
dx IS_A real_constant; |
| 36 |
dx :== 1.0 / (M - 1); |
| 37 |
|
| 38 |
coeff IS_A real_constant; |
| 39 |
coeff :== 1/dx/dx; |
| 40 |
|
| 41 |
u[0..M-1][0..M-1] IS_A solver_var; |
| 42 |
udot[0..M-1][0..M-1] IS_A solver_var; |
| 43 |
|
| 44 |
t IS_A time; |
| 45 |
|
| 46 |
(* central difference equations for internal points *) |
| 47 |
FOR i IN [1..M-2] CREATE |
| 48 |
FOR j IN [1..M-2] CREATE |
| 49 |
udoteq[i][j]: udot[i][j] = coeff * ( |
| 50 |
(* u_yy *) u[i-1][j] + u[i+1][j] - 2*u[i][j] |
| 51 |
+ (* u_xx *) u[i][j-1] + u[i][j+1] - 2*u[i][j] |
| 52 |
); |
| 53 |
END FOR; |
| 54 |
END FOR; |
| 55 |
|
| 56 |
(* boundary conditions: constrained u = 0 *) |
| 57 |
|
| 58 |
(* top and bottom *) |
| 59 |
FOR i IN [0,M-1] CREATE |
| 60 |
FOR j IN[0..M-1] CREATE |
| 61 |
topbot[i][j]: u[i][j] = 0; |
| 62 |
(* udot1[i][j]: udot[i][j] = 0; *) |
| 63 |
END FOR; |
| 64 |
END FOR; |
| 65 |
|
| 66 |
(* sides *) |
| 67 |
FOR i IN [1..M-2] CREATE |
| 68 |
FOR j IN[0,M-1] CREATE |
| 69 |
sides[i][j]: u[i][j] = 0; |
| 70 |
(* udot2[i][j]: udot[i][j] = 0; *) |
| 71 |
END FOR; |
| 72 |
END FOR; |
| 73 |
|
| 74 |
METHODS |
| 75 |
METHOD on_load; |
| 76 |
RUN specify; |
| 77 |
RUN values; |
| 78 |
RUN ode_init; |
| 79 |
END on_load; |
| 80 |
|
| 81 |
METHOD specify; |
| 82 |
FIX udot[1..M-2][1..M-2]; |
| 83 |
(* |
| 84 |
FIX udot[0,M-1][0..M-1]; |
| 85 |
FIX udot[1..M-2][0,M-1]; |
| 86 |
*) |
| 87 |
END specify; |
| 88 |
|
| 89 |
|
| 90 |
METHOD values; |
| 91 |
t := 0 {s}; |
| 92 |
FOR i IN [0..M-1] DO |
| 93 |
FOR j IN [0..M-1] DO |
| 94 |
u[i][j] := 16 * (dx*j) * (1 - dx*j) * dx*i * (1 - dx*i); |
| 95 |
END FOR; |
| 96 |
END FOR; |
| 97 |
FOR i IN [1..M-2] DO |
| 98 |
FOR j IN [1..M-2] DO |
| 99 |
u[i][j] := 0; |
| 100 |
END FOR; |
| 101 |
END FOR; |
| 102 |
END values; |
| 103 |
|
| 104 |
METHOD ode_init; |
| 105 |
FREE udot[1..M-2][1..M-2]; (* middle *) |
| 106 |
FIX udot[0,M-1][0..M-1]; (* top, bottom *) |
| 107 |
FIX udot[1..M-2][0,M-1]; (* sides *) |
| 108 |
|
| 109 |
FOR i IN [0..M-1] DO |
| 110 |
FOR j IN[0..M-1] DO |
| 111 |
u[i][j].ode_id := i*M + j; |
| 112 |
udot[i][j].ode_id := i*M + j; |
| 113 |
|
| 114 |
u[i][j].ode_type := 1; |
| 115 |
udot[i][j].ode_type := 2; |
| 116 |
END FOR; |
| 117 |
END FOR; |
| 118 |
|
| 119 |
t.obs_id := 1; |
| 120 |
udot[M/2][M/2].obs_id := 2; |
| 121 |
|
| 122 |
t.ode_type := -1; |
| 123 |
|
| 124 |
END ode_init; |
| 125 |
|
| 126 |
END idakryx; |
Parent Directory
| 
Revision Log