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