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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 968 - (hide annotations) (download) (as text)
Mon Dec 18 05:49:00 2006 UTC (13 years, 10 months ago) by johnpye
File MIME type: text/x-ascend
File size: 3014 byte(s)
Added SCons tests to check SIGINT and to replace ascresetneeded (need replacement for this in Autoconf as well).
Removed debugging from createinst.c
Typo (text) in evaluate.c
Commented out redundant code in importhandler.c
Added signal handling in ExecuteCASGN.
Added missing ospath_free in ModuleSearchPath.
Exported InitSymbolTable, DestroySymbolTable in symtab (dubious)
Moved FPRESET macro out of ascConfig.h and into ascSignal.h
Added Asc_SignalHandler{Push,Pop}Default.
Added ASC_RESETNEEDED and HAVE_C99FPE macros in config.h.in.
Found the bug causing the SIGFPE in idakryx.a4c (raises a question about int/float division in modelling, I think)
Added system_destroy call in Simulation::~Simulation (dubious).
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 956 M :== 4;
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     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