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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 968 - (show 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 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.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];
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