% fdm.s % % Purpose: Finite Difference Schemes % Copyright: R.A. van Engelen, Leiden University, 1997include 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.