/[ascend]/trunk/models/johnpye/idakryx.a4c
ViewVC logotype

Contents of /trunk/models/johnpye/idakryx.a4c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 956 - (show annotations) (download) (as text)
Sat Dec 9 15:38:05 2006 UTC (13 years, 10 months ago) by johnpye
File MIME type: text/x-ascend
File size: 3011 byte(s)
Removed references to 'ASC_BIG_BUGMAIL' and 'ASC_MILD_BUGMAIL'
Disabled stream redirection functions (they were making testing difficult)
Re-wrote the main test-runner to permit specific tests to be run from the commandline (ongoing)
Disabled readln tests that required stream redirection.
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 :== 4;
34
35 dx IS_A real_constant;
36 dx :== 1 / (M - 1);
37
38 coeff IS_A real_constant;
39 coeff :== 1/dx^2;
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];
106 FREE udot[0,M-1][0..M-1];
107 FREE udot[1..M-2][0,M-1];
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.ode_type := -1;
120
121 END ode_init;
122
123 END idakryx;

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22