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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 953 - (hide annotations) (download) (as text)
Thu Dec 7 14:47:15 2006 UTC (13 years, 10 months ago) by johnpye
File MIME type: text/x-ascend
File size: 2929 byte(s)
Added test for C99 FPE handling
Fixing mess-up of ChildByChar in arrayinst.h header.
Added 'safeeval' config option to IDA.
Changed 'SigHandler' to 'SigHandlerFn *' in line with other function pointer datatypes being used in ASCEND.
Moved processVarStatus *after* 'Failed integrator' exception (ongoing issue).
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     M :== 10;
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     END values;
98    
99     METHOD ode_init;
100     FREE udot[1..M-2][1..M-2];
101     FREE udot[0,M-1][0..M-1];
102     FREE udot[1..M-2][0,M-1];
103    
104     FOR i IN [0..M-1] DO
105     FOR j IN[0..M-1] DO
106     u[i][j].ode_id := i*M + j;
107     udot[i][j].ode_id := i*M + j;
108    
109     u[i][j].ode_type := 1;
110     udot[i][j].ode_type := 2;
111     END FOR;
112     END FOR;
113    
114     t.ode_type := -1;
115    
116     END ode_init;
117    
118     END idakryx;

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