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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 990 - (hide annotations) (download) (as text)
Thu Dec 21 07:12:58 2006 UTC (13 years, 10 months ago) by johnpye
File MIME type: text/x-ascend
File size: 3101 byte(s)
Added TestIDA.testkryxDENSE which seems to work OK.
1 johnpye 953 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 johnpye 990 M :== 5;
34 johnpye 953
35     dx IS_A real_constant;
36 johnpye 968 dx :== 1.0 / (M - 1);
37 johnpye 953
38     coeff IS_A real_constant;
39 johnpye 968 coeff :== 1/dx/dx;
40 johnpye 953
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 johnpye 956 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 johnpye 953 END values;
103    
104     METHOD ode_init;
105 johnpye 985 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 johnpye 953
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 johnpye 990 t.obs_id := 1;
120     udot[M/2][M/2].obs_id := 2;
121    
122 johnpye 953 t.ode_type := -1;
123    
124     END ode_init;
125    
126     END idakryx;

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