models/johnpye/advection.a4c

(*
	A model of the advection equation, "the simplest imaginable PDE,...
	the great difficulty of solving it numerically causes considerable
	embararassment to those attempting to solve complex PDE systems."

	Carver and Hinds, 'The method of lines and the advective equation'
	/Simulation/, 1978.
	http://sim.sagepub.com/cgi/content/abstract/31/2/59

	(incidentally, in that paper, watch out for the error in equation (13) in 
	the sign of the u_{i-2} term)

	Here we model the linear advection equation as given by Carver and Hinds.
	We will investigate the affect of different spatial discretisations when
	integrating with the multi-step integrators (IDA and LSODE) currently 
	provided with ASCEND.

	by John Pye, Feb 2007.
*)
REQUIRE "ivpsystem.a4l";
REQUIRE "atoms.a4l";

MODEL advection_common;

	(* discretisation *)
	L IS_A real_constant;
	L :== 10 {m};
	n IS_A integer_constant;
	n :== 200;
	dx IS_A real_constant;
	dx :== L / n;

	(* node sets *)
	nodes, blip IS_A set OF integer_constant;
	nodes :== [1..n];
	blip :== [10..29];

	(* advection speed *)
	v IS_A speed;

	(* the advecting quantity *)
	u[nodes] IS_A mass_density;
	du_dt[nodes] IS_A density_rate;
	t IS_A time;

	x[nodes] IS_A real_constant;
	FOR i IN nodes CREATE
		x[i] :== i * dx;
	END FOR;

METHODS
	METHOD default;
		FOR i IN nodes DO
			u[i] := 1 {kg/m^3};
			du_dt[i] := 0 {kg/m^3/s};
		END FOR;
	END default;

	METHOD specify;
		FIX du_dt[1];
		FIX u[1];
		FIX v;
	END specify;
	
	METHOD values;
		t := 0 {s};
		v := 1 {m/s};
		du_dt[1] := 0 {kg/m^3/s};
		u[1] := 1 {kg/m^3};
		u[blip] := 10 {kg/m^3};
	END values;

	METHOD bound_self;
		u[nodes].lower_bound := -1e5 {kg/m^3};
		u[nodes].upper_bound := 1e5 {kg/m^3};
	END bound_self;

	METHOD on_load;
		RUN default_self;
		RUN bound_self;
		RUN ode_init;
		RUN reset;
		RUN values;
	END on_load;

	METHOD reset_steady;
		RUN ClearAll;
		RUN specify;
		FIX du_dt[butfirst];
	END reset_steady;

	METHOD ode_init;
		FOR i IN nodes DO
			u[i].ode_id := i;
			du_dt[i].ode_id := i;
			u[i].ode_type := 1;
			du_dt[i].ode_type := 2;
		END FOR;
		t.ode_type := -1;
		
		FOR i IN [5..120] DO
			u[i].obs_id := 2*i;
		END FOR;

	END ode_init;

END advection_common;

MODEL advection_upwind2 REFINES advection_common;

	upwind2 IS_A set OF integer_constant;
	upwind2 :== nodes - [1];

	(* use two-point upwind approximation to du_dx *)
	FOR i IN upwind2 CREATE
		du_dt[i] + v * (u[i] - u[i - 1]) / dx = 0;
	END FOR;

END advection_upwind2;

MODEL advection_central3 REFINES advection_common;

	central3 IS_A set OF integer_constant;
	central3 :== nodes - [1,n];

	(* use two-point upwind approximation for du_dt[2] and du_dt[n] *)
	FOR i IN [n] CREATE
		du_dt[i] + v * (u[i] - u[i - 1]) / dx = 0;
	END FOR;

	(* use central difference 3-point approximation for remaining points *)
	FOR i IN central3 CREATE
		du_dt[i] + v * (u[i+1] - u[i-1]) / 2. / dx = 0;
	END FOR;

END advection_central3;

MODEL advection_upwind3 REFINES advection_common;

	upwind3 IS_A set OF integer_constant;
	upwind3 :== nodes - [1,2];

	(* use two-point upwind approximation for du_dt[2] *)
	FOR i IN [2] CREATE
		du_dt[i] + v * (u[i] - u[i - 1]) / dx = 0;
	END FOR;

	(* use three-point upwind approximation for remaining points *)
	FOR i IN upwind3 CREATE
		du_dt[i] + v * (3 * u[i] - 4 * u[i-1] + u[i-2]) / 2. / dx = 0;
	END FOR;

END advection_upwind3;

MODEL advection_upwind4 REFINES advection_common;

	upwind4 IS_A set OF integer_constant;
	upwind4 :== nodes - [1,2,n];

	(* use two-point upwind approximation for du_dt[2] and du_dt[n]*)
	FOR i IN [2,n] CREATE
		du_dt[i] + v * (u[i] - u[i - 1]) / dx = 0;
	END FOR;

	(* use four-point upwind approximation for remaining points *)
	FOR i IN upwind4 CREATE
		du_dt[i] + v * (-u[i+1] + 6.*u[i] - 3.*u[i-1] - 2.*u[i-2]) / 6. / dx = 0;
	END FOR;

END advection_upwind4;
1	(*
2	A model of the advection equation, "the simplest imaginable PDE,...
3	the great difficulty of solving it numerically causes considerable
4	embararassment to those attempting to solve complex PDE systems."
5
6	Carver and Hinds, 'The method of lines and the advective equation'
7	/Simulation/, 1978.
8	http://sim.sagepub.com/cgi/content/abstract/31/2/59
9
10	(incidentally, in that paper, watch out for the error in equation (13) in
11	the sign of the u_{i-2} term)
12
13	Here we model the linear advection equation as given by Carver and Hinds.
14	We will investigate the affect of different spatial discretisations when
15	integrating with the multi-step integrators (IDA and LSODE) currently
16	provided with ASCEND.
17
18	by John Pye, Feb 2007.
19	*)
20	REQUIRE "ivpsystem.a4l";
21	REQUIRE "atoms.a4l";
22
23	MODEL advection_common;
24
25	(* discretisation *)
26	L IS_A real_constant;
27	L :== 10 {m};
28	n IS_A integer_constant;
29	n :== 200;
30	dx IS_A real_constant;
31	dx :== L / n;
32
33	(* node sets *)
34	nodes, blip IS_A set OF integer_constant;
35	nodes :== [1..n];
36	blip :== [10..29];
37
38	(* advection speed *)
39	v IS_A speed;
40
41	(* the advecting quantity *)
42	u[nodes] IS_A mass_density;
43	du_dt[nodes] IS_A density_rate;
44	t IS_A time;
45
46	x[nodes] IS_A real_constant;
47	FOR i IN nodes CREATE
48	x[i] :== i * dx;
49	END FOR;
50
51	METHODS
52	METHOD default;
53	FOR i IN nodes DO
54	u[i] := 1 {kg/m^3};
55	du_dt[i] := 0 {kg/m^3/s};
56	END FOR;
57	END default;
58
59	METHOD specify;
60	FIX du_dt[1];
61	FIX u[1];
62	FIX v;
63	END specify;
64
65	METHOD values;
66	t := 0 {s};
67	v := 1 {m/s};
68	du_dt[1] := 0 {kg/m^3/s};
69	u[1] := 1 {kg/m^3};
70	u[blip] := 10 {kg/m^3};
71	END values;
72
73	METHOD bound_self;
74	u[nodes].lower_bound := -1e5 {kg/m^3};
75	u[nodes].upper_bound := 1e5 {kg/m^3};
76	END bound_self;
77
78	METHOD on_load;
79	RUN default_self;
80	RUN bound_self;
81	RUN ode_init;
82	RUN reset;
83	RUN values;
84	END on_load;
85
86	METHOD reset_steady;
87	RUN ClearAll;
88	RUN specify;
89	FIX du_dt[butfirst];
90	END reset_steady;
91
92	METHOD ode_init;
93	FOR i IN nodes DO
94	u[i].ode_id := i;
95	du_dt[i].ode_id := i;
96	u[i].ode_type := 1;
97	du_dt[i].ode_type := 2;
98	END FOR;
99	t.ode_type := -1;
100
101	FOR i IN [5..120] DO
102	u[i].obs_id := 2*i;
103	END FOR;
104
105	END ode_init;
106
107	END advection_common;
108
109	MODEL advection_upwind2 REFINES advection_common;
110
111	upwind2 IS_A set OF integer_constant;
112	upwind2 :== nodes - [1];
113
114	(* use two-point upwind approximation to du_dx *)
115	FOR i IN upwind2 CREATE
116	du_dt[i] + v * (u[i] - u[i - 1]) / dx = 0;
117	END FOR;
118
119	END advection_upwind2;
120
121	MODEL advection_central3 REFINES advection_common;
122
123	central3 IS_A set OF integer_constant;
124	central3 :== nodes - [1,n];
125
126	(* use two-point upwind approximation for du_dt[2] and du_dt[n] *)
127	FOR i IN [n] CREATE
128	du_dt[i] + v * (u[i] - u[i - 1]) / dx = 0;
129	END FOR;
130
131	(* use central difference 3-point approximation for remaining points *)
132	FOR i IN central3 CREATE
133	du_dt[i] + v * (u[i+1] - u[i-1]) / 2. / dx = 0;
134	END FOR;
135
136	END advection_central3;
137
138	MODEL advection_upwind3 REFINES advection_common;
139
140	upwind3 IS_A set OF integer_constant;
141	upwind3 :== nodes - [1,2];
142
143	(* use two-point upwind approximation for du_dt[2] *)
144	FOR i IN [2] CREATE
145	du_dt[i] + v * (u[i] - u[i - 1]) / dx = 0;
146	END FOR;
147
148	(* use three-point upwind approximation for remaining points *)
149	FOR i IN upwind3 CREATE
150	du_dt[i] + v * (3 * u[i] - 4 * u[i-1] + u[i-2]) / 2. / dx = 0;
151	END FOR;
152
153	END advection_upwind3;
154
155	MODEL advection_upwind4 REFINES advection_common;
156
157	upwind4 IS_A set OF integer_constant;
158	upwind4 :== nodes - [1,2,n];
159
160	(* use two-point upwind approximation for du_dt[2] and du_dt[n]*)
161	FOR i IN [2,n] CREATE
162	du_dt[i] + v * (u[i] - u[i - 1]) / dx = 0;
163	END FOR;
164
165	(* use four-point upwind approximation for remaining points *)
166	FOR i IN upwind4 CREATE
167	du_dt[i] + v * (-u[i+1] + 6.u[i] - 3.u[i-1] - 2.*u[i-2]) / 6. / dx = 0;
168	END FOR;
169
170	END advection_upwind4;