% 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.