%	fdm.s
%
%	Purpose:	Finite Difference Schemes
%	Copyright:	R.A. van Engelen, Leiden University, 1997

include equation.

space (x(i), y(j), z(k)) time t(n).

(X :: T) ~ (_G2512 :: string) :: T := X.

x :: coordinate ~ "m".
y :: coordinate ~ "m".
z :: coordinate ~ "m".
t :: coordinate ~ "s".

grid is grid_type.
half is grid_type.

x(i) :: integer ~ "m" field x(grid) := i.
y(j) :: integer ~ "m" field y(grid) := j.
z(k) :: integer ~ "m" field z(grid) := k.

prefix delta has latex_name("\Delta").

delta (_G2512 :: coordinate ~ Unit) :: real ~ Unit.

index(_G2509) is converter_of(type) and elemental.

index(x) :: index := i.
index(y) :: index := j.
index(z) :: index := k.

point(_G2509) is converter_of(type) and elemental.

point((X :: coordinate)) :: integer := index(X).

index_reference(_G2509) is converter_of(type) and elemental.

(index_reference((X = Y :: reference(coordinate))) ::
reference(index) := index(X) = Y)
.
index_domain(_G2509) is converter_of(type) and elemental. (index_domain((X = R :: domain(coordinate))) ::
domain(index) := index(X) = R)
.
((_G2512 :: T ~ U) @ (_G2526 :: reference(index) ~ "1" field F) ::
T ~ U field F)
.
E @ R latex precedence(none); "\left[", E, "\right]_", R. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Grid-location protection % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% protect((X :: T ~ U field x), x) :: T ~ U field x := X. protect((X :: T ~ U field y), y) :: T ~ U field y := X. protect((X :: T ~ U field z), z) :: T ~ U field z := X. (protect(X, C) latex
precedence(none); "\overbrace", X, "^{{\rm protect}(", C, ")}")
.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Relative (half-) grid-point references % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% shift((_G2512 :: T ~ Unit1), (_G2518 :: coordinate ~ Unit2), (_G2524 :: integer ~ Unit2)) :: T ~ Unit1. shift((_G2512 :: T ~ Unit1), (_G2518 :: coordinate ~ Unit2), (_G2524 :: real ~ Unit2)) :: T ~ Unit1. (shift(_G2509, _G2510, _G2511) is
(stencil_op and
linear and
abstraction_of(grid, choices = ( shift_g
, shift_h
, shift_c
))
))
.
shift_g(_G2509, _G2510, _G2511) is stencil_op and linear and self_commuting. (shift_h(_G2509, _G2510, _G2511) is
(stencil_op and linear and self_commuting can
commute_over(shift_g arity 3)
))
.
(shift_c(_G2509, _G2510, _G2511) is
(stencil_op and linear and self_commuting can
commute_over(shift_g arity 3, shift_h arity 3)
))
.
(shift_g((E :: T field x(grid)), x, (- N / 2 :: real)) ::
T field x(half) := E @ (index(x) = point(x) - floor(N / 2)))
.
(shift_g((E :: T field x(grid)), x, (N / 2 :: real)) ::
T field x(half) := E @ (index(x) = point(x) + floor(N / 2)))
.
(shift_h((E :: T field x(half)), x, (- N / 2 :: real)) ::
T field x(grid) := E @ (index(x) = point(x) - floor((N + 1) / 2)))
.
(shift_h((E :: T field x(half)), x, (N / 2 :: real)) ::
T field x(grid) := E @ (index(x) = point(x) + floor((N + 1) / 2)))
.
(shift_c((E :: T field x(G)), x, (N :: integer)) ::
T field x(G) := E @ (index(x) = point(x) + floor(N / 2)))
.
(shift_g((E :: T field y(grid)), y, (- N / 2 :: real)) ::
T field y(half) := E @ (index(y) = point(y) - floor(N / 2)))
.
(shift_g((E :: T field y(grid)), y, (N / 2 :: real)) ::
T field y(half) := E @ (index(y) = point(y) + floor(N / 2)))
.
(shift_h((E :: T field y(half)), y, (- N / 2 :: real)) ::
T field y(grid) := E @ (index(y) = point(y) - floor((N + 1) / 2)))
.
(shift_h((E :: T field y(half)), y, (N / 2 :: real)) ::
T field y(grid) := E @ (index(y) = point(y) + floor((N + 1) / 2)))
.
(shift_c((E :: T field y(G)), y, (N :: integer)) ::
T field y(G) := E @ (index(y) = point(y) + floor(N / 2)))
.
(shift_g((E :: T field z(grid)), z, (- N / 2 :: real)) ::
T field z(half) := E @ (index(z) = point(z) - floor(N / 2)))
.
(shift_g((E :: T field z(grid)), z, (N / 2 :: real)) ::
T field z(half) := E @ (index(z) = point(z) + floor(N / 2)))
.
(shift_h((E :: T field z(half)), z, (- N / 2 :: real)) ::
T field z(grid) := E @ (index(z) = point(z) - floor((N + 1) / 2)))
.
(shift_h((E :: T field z(half)), z, (N / 2 :: real)) ::
T field z(grid) := E @ (index(z) = point(z) + floor((N + 1) / 2)))
.
(shift_c((E :: T field z(G)), z, (N :: integer)) ::
T field z(G) := E @ (index(z) = point(z) + floor(N / 2)))
.
% NOTE: the index(_) and point(_) functions above are necessary, because type checking failes due to the variable % type T for the first argument E of shift. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Finite difference stencils % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% (df((_G2512 :: real ~ Unit1), (_G2518 :: coordinate ~ Unit2)) ::
real ~ (Unit1 / Unit2))
.
(df(_G2509, _G2510) is
(stencil_op and
differentiation_op and
abstraction_of(grid, choices = (df_g, df_h, df_c))
))
.
df_g arity 2 is stencil_op and differentiation_op. df_h arity 2 is stencil_op and differentiation_op. df_c arity 2 is stencil_op and differentiation_op. (df_g((E :: real ~ Unit1), (X :: coordinate ~ Unit2)) ::
real ~ (Unit1 / Unit2) := (E @ (X = X + 1) - E) / delta X)
.
(df_h((E :: real ~ Unit1), (X :: coordinate ~ Unit2)) ::
real ~ (Unit1 / Unit2) := (E - E @ (X = X - 1)) / delta X)
.
(df_c((E :: real ~ Unit1), (X :: coordinate ~ Unit2)) ::
(real ~ (Unit1 / Unit2) :=
(E @ (X = X + 1) - E @ (X = X - 1)) / (2 * delta X)
))
.
df_g((_G2513 :: field x(grid)), x) :: field x(half). df_h((_G2513 :: field x(half)), x) :: field x(grid). df_c((_G2513 :: field x(grid)), x) :: field x(grid). df_g((_G2513 :: field y(grid)), y) :: field y(half). df_h((_G2513 :: field y(half)), y) :: field y(grid). df_c((_G2513 :: field y(grid)), y) :: field y(grid). df_g((_G2513 :: field z(grid)), z) :: field z(half). df_h((_G2513 :: field z(half)), z) :: field z(grid). df_c((_G2513 :: field z(grid)), z) :: field z(grid). df_g((_G2513 :: field t(grid)), t) :: field t(half). df_h((_G2513 :: field t(half)), t) :: field t(grid). df_c((_G2513 :: field t(grid)), t) :: field t(grid). df_g(E, X) latex precedence(300); "\delta^+_", X, E. df_h(E, X) latex precedence(300); "\delta^-_", X, E. df_c(E, X) latex precedence(300); "\delta_", X, E. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Integration % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% (sum((_G2515 :: real ~ U field F), (_G2526 :: domain(index) ~ "1")) ::
real ~ U field F)
.
(int(_G2509, _G2510) is
(stencil_op and
integration_op and
abstraction_of(grid, choices = (mid_g, mid_h, trap))
))
.
(int((_G2512 :: real ~ Unit1), (_G2520 :: domain(coordinate) ~ Unit2)) ::
real ~ (Unit1 * Unit2))
.
mid_g(_G2509, _G2510) is stencil_op and integration_op. (mid_g((E :: real), (X = L .. U :: domain(coordinate))) ::
real := sum(delta X * E, X = L + 1 .. U))
.
mid_g((_G2513 :: field x(grid)), x = _G2516 .. _G2517) :: field x(half). mid_g((_G2513 :: field y(grid)), y = _G2516 .. _G2517) :: field y(half). mid_g((_G2513 :: field z(grid)), z = _G2516 .. _G2517) :: field z(half). mid_h(_G2509, _G2510) is stencil_op and integration_op. (mid_h((E :: real), (X = L .. U :: domain(coordinate))) ::
real := sum(delta X * E, X = L .. U - 1))
.
mid_h((_G2513 :: field x(half)), x = _G2516 .. _G2517) :: field x(grid). mid_h((_G2513 :: field y(half)), y = _G2516 .. _G2517) :: field y(grid). mid_h((_G2513 :: field z(half)), z = _G2516 .. _G2517) :: field z(grid). trap(_G2509, _G2510) is stencil_op and integration_op. (trap((E :: real field F), (X = L .. U :: domain(coordinate))) ::
(real field F :=
( sum(delta X * E, X = L + 1 .. U - 1)
+ 1 / 2 * delta X * (E @ (X = L) + E @ (X = U))
)
))
.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Location % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% firstloc((_G2512 :: boolean field F), (_G2523 :: domain(index) ~ "1"), (_G2529 :: integer field F)) :: integer field F. lastloc((_G2512 :: boolean field F), (_G2523 :: domain(index) ~ "1"), (_G2529 :: integer field F)) :: integer field F.