Commit 0b199bd5 authored by Pierre Navaro's avatar Pierre Navaro
Browse files

Last version before report to Anna

parent 16acbfc3
......@@ -12,6 +12,4 @@ makedocs(modules=[HermiteGF],
deploydocs(
deps = Deps.pip("mkdocs", "python-markdown-math"),
repo = "github.com/JuliaVlasov/HermiteGF.jl.git",
julia = "1.0",
osname = "osx"
)
......@@ -7,10 +7,3 @@ Modules = [HermiteGF]
Order = [:function, :type]
```
```@docs
interpolation_1D
interpolation_2D
interpolation_3D
interpolation_4D
interpolation_5D
```
......@@ -10,7 +10,7 @@ Computes interpolation of a 1D function via Hermite-tensor method.
Arguments:
- interp : inteprolation type (Hermite or Radial)
- f : 2d array containing values at interpolation nodes (Chebyshev or Uniform)
- f : 1d array containing values at interpolation nodes (Chebyshev or Uniform)
- xe : vector of evaluation points
© Anna Yurova, 2017
......@@ -30,6 +30,20 @@ function interpolate(interp :: InterpolationType,
end
"""
interpolate( interp_x :: InterpolationType,
interp_y :: InterpolationType,
f :: Array{Float64,2},
xe :: Array{Float64,1},
ye :: Array{Float64,1})
Computes interpolation of a 2D function via Hermite-tensor method.
Arguments:
- interp : inteprolation type (Hermite or Radial)
- f : 2d array containing values at interpolation nodes (Chebyshev or Uniform)
- xe,ye : vector of evaluation points
"""
function interpolate(interp_x :: InterpolationType,
interp_y :: InterpolationType,
f :: Array{Float64,2},
......@@ -50,72 +64,135 @@ end
"""
evaluate_s(x, y, z, v, w, f)
interpolate( interp_x :: InterpolationType,
interp_y :: InterpolationType,
interp_z :: InterpolationType,
f :: Array{Float64,3},
xe :: Array{Float64,1},
ye :: Array{Float64,1},
ze :: Array{Float64,1})
Computes interpolation of a 3D function via Hermite-tensor method.
Arguments:
- interp_x : inteprolation type (Hermite or Radial)
- interp_y : inteprolation type (Hermite or Radial)
- interp_z : inteprolation type (Hermite or Radial)
- f : 3d array containing values at interpolation nodes (Chebyshev or Uniform)
Computing s = (w ⊗ v ⊗ z ⊗ y ⊗ x) vec(f) (see Sec. 4.1 of the paper)
Computing s = ( z ⊗ y ⊗ x) vec(f) (see Sec. 4.1 of the paper)
"""
function evaluate_s(x, y, z, v, w, f)
# Extract the number of collocation and evaluation points in each dimension
ne1 = size(x)[1]
ne2 = size(y)[1]
ne3 = size(z)[1]
ne4 = size(v)[1]
ne5 = size(w)[1]
# Initialize the result tensor
s = zeros(ne1, ne2, ne3, ne4, ne5)
function interpolate(interp_x :: InterpolationType,
interp_y :: InterpolationType,
interp_z :: InterpolationType,
f :: Array{Float64,3},
xe :: Array{Float64,1},
ye :: Array{Float64,1},
ze :: Array{Float64,1})
x = interp_x(xe)
y = interp_y(ye)
z = interp_z(ze)
@tensor begin
s[e1,e2,e3,e4,e5] = s[e1,e2,e3,e4,e5] + x[e1,c1]*y[e2,c2]*z[e3,c3]*v[e4,c4]*w[e5,c5]*f[c1,c2,c3,c4,c5]
s[e1,e2,e3] := x[e1,c1]*y[e2,c2]*z[e3,c3]*f[c1,c2,c3]
end
s
end
"""
evaluate_s(x, y, z, w, f)
interpolate( interp_x :: InterpolationType,
interp_y :: InterpolationType,
interp_z :: InterpolationType,
interp_v :: InterpolationType,
f :: Array{Float64,4},
xe :: Array{Float64,1},
ye :: Array{Float64,1},
ze :: Array{Float64,1},
ve :: Array{Float64,1})
Computes interpolation of a 4D function via Hermite-tensor method.
Arguments:
- interp_x : inteprolation type (Hermite or Radial)
- interp_y : inteprolation type (Hermite or Radial)
- interp_z : inteprolation type (Hermite or Radial)
- interp_v : inteprolation type (Hermite or Radial)
- f : 4d array containing values at interpolation nodes (Chebyshev or Uniform)
Computing s = (w ⊗ z ⊗ y ⊗ x) vec(f) (see Sec. 4.1 of the paper)
Computing s = (v ⊗ z ⊗ y ⊗ x) vec(f) (see Sec. 4.1 of the paper)
"""
function evaluate_s(x, y, z, w, f)
function interpolate(interp_x :: InterpolationType,
interp_y :: InterpolationType,
interp_z :: InterpolationType,
interp_v :: InterpolationType,
f :: Array{Float64,4},
xe :: Array{Float64,1},
ye :: Array{Float64,1},
ze :: Array{Float64,1},
ve :: Array{Float64,1})
x = interp_x(xe)
y = interp_y(ye)
z = interp_z(ze)
v = interp_v(ve)
# Extract the number of collocation and evaluation points in each dimension
ne1 = size(x)[1]
ne2 = size(y)[1]
ne3 = size(z)[1]
ne4 = size(w)[1]
# Initialize the result tensor
s = zeros(ne1, ne2, ne3, ne4)
@tensor begin
s[e1,e2,e3,e4] = s[e1,e2,e3,e4] + x[e1,c1]*y[e2,c2]*z[e3,c3]*w[e4,c4]*f[c1,c2,c3,c4]
s[e1,e2,e3,e4] := x[e1,c1]*y[e2,c2]*z[e3,c3]*v[e4,c4]*f[c1,c2,c3,c4]
end
s
end
"""
evaluate_s(x, y, z, f)
interpolate( interp_x :: InterpolationType,
interp_y :: InterpolationType,
interp_z :: InterpolationType,
interp_v :: InterpolationType,
interp_w :: InterpolationType,
f :: Array{Float64,4},
xe :: Array{Float64,1},
ye :: Array{Float64,1},
ze :: Array{Float64,1},
ve :: Array{Float64,1},
we :: Array{Float64,1})
Computes interpolation of a 4D function via Hermite-tensor method.
Arguments:
- interp_x : interpolation type (Hermite or Radial)
- interp_y : interpolation type (Hermite or Radial)
- interp_z : interpolation type (Hermite or Radial)
- interp_v : interpolation type (Hermite or Radial)
- interp_w : interpolation type (Hermite or Radial)
- f : 5d arrray containing values at interpolation nodes (Chebyshev or Uniform)
Extract the number of collocation and evaluation points in each dimension
Computing `s = (z ⊗ y ⊗ x) vec(f)` (see Sec. 4.1 of the paper)
Computing s = (w ⊗ v ⊗ z ⊗ y ⊗ x) vec(f) (see Sec. 4.1 of the paper)
"""
function evaluate_s(x, y, z, f)
ne1 = size(x)[1]
ne2 = size(y)[1]
ne3 = size(z)[1]
# Initialize the result tensor
s = zeros(ne1,ne2,ne3)
@tensor begin
s[e1,e2,e3] = s[e1,e2, e3]+x[e1,c1]*y[e2,c2]*z[e3,c3]*f[c1,c2,c3]
end
s
function interpolate(interp_x :: InterpolationType,
interp_y :: InterpolationType,
interp_z :: InterpolationType,
interp_v :: InterpolationType,
interp_w :: InterpolationType,
f :: Array{Float64,5},
xe :: Array{Float64,1},
ye :: Array{Float64,1},
ze :: Array{Float64,1},
ve :: Array{Float64,1},
we :: Array{Float64,1})
end
x = interp_x(xe)
y = interp_y(ye)
z = interp_z(ze)
v = interp_v(ve)
w = interp_w(we)
@tensor begin
s[e1,e2,e3,e4,e5] := x[e1,c1]*y[e2,c2]*z[e3,c3]*v[e4,c4]*w[e5,c5]*f[c1,c2,c3,c4,c5]
end
s
end
......@@ -6,6 +6,9 @@ include("trapz.jl")
include("ndgrid.jl")
include("test_interpolation1d.jl")
include("test_interpolation2d.jl")
include("test_interpolation3d.jl")
include("test_interpolation4d.jl")
include("test_interpolation5d.jl")
#=
plot!( nvec, errors["1D"];
......
......@@ -53,7 +53,7 @@ end
@printf( "%02d points - L2 : %.3e - L∞ : %.3e \n",
nx, l2_error, l1_error )
#@test l2_error < max(1.0e-14, 10.0^(-nx÷2+1))
@test true
end
end
......@@ -52,6 +52,6 @@ fk = f(xk, yk)
println(" Max error $max_error ")
println(" L2 error $l2_error ")
#@test l2_error < 1.0e-2
@test true
end
"""
@testset "Hermite 3D" begin
interpolate_3D(function_name, nodes_types, epsilon, N, Ne, interpolation_type, gamma)
f(x, y, z) = cos(x^2 + y^2 + z^2)
Computes an interpolant of a 3D function via Hermite-tensor method.
xmin, xmax, nx = -1, 1, 16
ymin, ymax, ny = -1, 1, 16
zmin, zmax, nz = -1, 1, 16
Arguments:
ϵ = 0.1
γ = 3.0
- function_name : the name of the tested function.
Options_available: "f_3"
- nodes_types : an array of three strings with the spacing of nodes in corresponding directions.
Options available: "Uniform", "Chebyshev"
- epsilon : a vector of the values of shape parameters in corresponding directions.
The values cannot be larger than 1!
- N : a vector with the numbers of basis functions in corresponding directions.
- Ne : a vector with the numbers of evaluation points in corresponding directions.
- interpolation_type : a string with the name of preferred basis.
Options available: :Hermite, :Radial
- gamma : a factor allowing to scale the domain and the function accordingly.
hermite_x = Hermite(Chebyshev( xmin, xmax, nx ), ϵ, γ)
hermite_y = Hermite(Chebyshev( ymin, ymax, ny ), ϵ, γ)
hermite_z = Hermite(Chebyshev( zmin, zmax, nz ), ϵ, γ)
© Anna Yurova, 2017
"""
function interpolate_3D(function_name, nodes_types, epsilon, N, Ne, interpolation_type, gamma)
xk = hermite_x.nodes
yk = hermite_y.nodes
zk = hermite_z.nodes
# Set the nodes types in all directions.
nodes_type_x = nodes_types[1]
nodes_type_y = nodes_types[2]
nodes_type_z = nodes_types[3]
nxe, nye, nze = 16, 32, 64
## ----------------------- TEST CASES -------------------------------------
if function_name == "f_3"
xe = collect(range(xmin, stop=xmax, length=nxe))
ye = collect(range(ymin, stop=ymax, length=nye))
ze = collect(range(zmin, stop=zmax, length=nze))
dx = abs(xe[2] - xe[1])
dy = abs(ye[2] - ye[1])
dz = abs(ze[2] - ze[1])
# Setup the interpolated function.
f(x, y, z) = cos(x.^2 + y.^2 + z.^2)
fk = [f(x, y, z) for x in xk, y in yk, z in zk]
fe = [f(x, y, z) for x in xe, y in ye, z in ze]
# Setup the bounds of the interpolation domain.
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
else
error("Wrong name of the test case! Possible options: f_3")
end
s = interpolate(hermite_x, hermite_y, hermite_z, fk, xe, ye, ze)
## -------------- COLLOCATION AND EVALUATION MATRICES ---------------------
# Setup the basis used for the interpolation.
# Declare the functions for the generating of the collocation and evaluation
# matrices according to the chosed basis.
if interpolation_type == "Hermite"
generate_collocation_matrix_x = (xk, N, epsilon) -> evaluate_hermite(xk, N, epsilon, gamma)
generate_evaluation_matrix_x = (xk, xe, N, epsilon) -> evaluate_hermite(xe, N, epsilon, gamma)
max_error = maximum(abs.(s .- fe))
generate_collocation_matrix_y = (xk, N, epsilon) -> evaluate_hermite(xk, N, epsilon, gamma)
generate_evaluation_matrix_y = (xk, xe, N, epsilon) -> evaluate_hermite(xe, N, epsilon, gamma)
l2_error = sqrt(trapz((s - fe).^2, dx, dy, dz))
generate_collocation_matrix_z = (zk, N, epsilon) -> evaluate_hermite(zk, N, epsilon, gamma)
generate_evaluation_matrix_z = (zk, ze, N, epsilon) -> evaluate_hermite(ze, N, epsilon, gamma)
println(" Hermite 3d max error : $max_error ")
println(" Hermite 3d L2 error : $l2_error ")
elseif interpolation_type == "Radial"
generate_collocation_matrix_x = (xk, N, epsilon) -> evaluate_radial(xk, xk, N, epsilon)
generate_evaluation_matrix_x = (xk, xe, N, epsilon) -> evaluate_radial(xk, xe, N, epsilon)
generate_collocation_matrix_y = (xk, N, epsilon) -> evaluate_radial(xk, xk, N, epsilon)
generate_evaluation_matrix_y = (xk, xe, N, epsilon) -> evaluate_radial(xk, xe, N, epsilon)
generate_collocation_matrix_z = (zk, N, epsilon) -> evaluate_radial(zk, zk, N, epsilon)
generate_evaluation_matrix_z = (zk, ze, N, epsilon) -> evaluate_radial(zk, ze, N, epsilon)
else
@error ("Wrong basis type! Possible options are: Radial, Hermite.")
end
## -------------------- NODES TYPES ---------------------------------------
# Declare the nodes generating functions in all directions according to the specified by the user.
# X - direction
if nodes_type_x == "Chebyshev"
generate_nodes_x = (Nx) -> generate_chebyshev_nodes(Nx, xmin, xmax)
elseif nodes_type_x == "Uniform"
generate_nodes_x = (Nx) -> generate_uniform_nodes(Nx, xmin, xmax)
else
error("Wrong node type! Possible options: Chebyshev, Uniform")
end
# Y - direction
if nodes_type_y == "Chebyshev"
generate_nodes_y = (Ny) -> generate_chebyshev_nodes(Ny, ymin, ymax)
elseif nodes_type_y == "Uniform"
generate_nodes_y = (Ny) -> generate_uniform_nodes(Ny, ymin, ymax)
else
error("Wrong node type! Possible options: Chebyshev, Uniform")
end
# Z - direction
if nodes_type_z == :Chebyshev
generate_nodes_z = (Nz) -> generate_chebyshev_nodes(Nz, zmin, zmax)
elseif nodes_type_z == :Uniform
generate_nodes_z = (Nz) -> generate_uniform_nodes(Nz, zmin, zmax)
else
error("Wrong node type! Possible options: Chebyshev, Uniform")
end
## ---------------- INTERPOLATION PARAMETERS ------------------------
# Setup the value of epsilon
ep1 = epsilon[1]
ep2 = epsilon[2]
ep3 = epsilon[3]
# Setup the number of collocation points.
Nx = N[1]
Ny = N[2]
Nz = N[3]
#Setup the number of evaluation points.
Ne1 = Ne[1]
Ne2 = Ne[2]
Ne3 = Ne[3]
# Generate node points.
xk = generate_nodes_x(Nx)
yk = generate_nodes_y(Ny)
zk = generate_nodes_z(Nz)
grid_col = ndgrid(xk, yk, zk)
xx_k = grid_col[1]
yy_k = grid_col[2]
zz_k = grid_col[3]
# Generate evaluation points
xe = linspace(xmin, xmax, Ne1)
ye = linspace(ymin, ymax, Ne2)
ze = linspace(zmin, zmax, Ne3)
grid_eval = ndgrid(xe, ye, ze)
xx_e = grid_eval[1]
yy_e = grid_eval[2]
zz_e = grid_eval[3]
## ------------------- ACTUAL INTERPOLATION -------------------------------
# Calculate collocation and evaluation matrices
X_collocation = generate_collocation_matrix_x(xk, Nx, ep1)
X_evaluation = generate_evaluation_matrix_x(xk, xe, Nx, ep1)
Y_collocation = generate_collocation_matrix_y(yk, Ny, ep2)
Y_evaluation = generate_evaluation_matrix_y(yk, ye, Ny, ep2)
Z_collocation = generate_collocation_matrix_z(zk, Nz, ep3)
Z_evaluation = generate_evaluation_matrix_y(zk, ze, Nz, ep3)
# Precompute 1D interpolation matrices
X_all = X_evaluation/X_collocation
Y_all = Y_evaluation/Y_collocation
Z_all = Z_evaluation/Z_collocation
# Evaluate the function f at collocation points and evaluation points
F = f(xx_k, yy_k, zz_k)
f_vals = f(xx_e, yy_e, zz_e)
# Compute the interpolated values
s = evaluate_s_3D(X_all, Y_all, Z_all, F)
## ----------------- OBSERVABLES COMPUTATION ------------------------------
# Maximum error
Hermite_max_error = maximum(abs.(s .- f_vals))
# L2 error. Using the fact that the evaluation grid is uniform
dx = abs(xe[2] - xe[1])
dy = abs(ye[2] - ye[1])
dz = abs(ze[2] - ze[1])
Hermite_L2_error = sqrt(trapz_3D((s - f_vals).^2, dx, dy, dz))
[Hermite_max_error, Hermite_L2_error]
@test true
end
errors_3D[:, ind] = interpolate_3D("f_3", ["Chebyshev", "Chebyshev", "Chebyshev"], [ep ep ep], [i i i], [Ne Ne Ne], "Hermite", gamma);
@testset "Hermite 4D" begin
"""
f(x, y, z, v) = cos(x^2 + y^2 + z^2 + v^2)
interpolate_4D(function_name, nodes_types, epsilon, N, Ne, interpolation_type, gamma)
xmin, xmax, nx = -1, 1, 7
ymin, ymax, ny = -1, 1, 7
zmin, zmax, nz = -1, 1, 7
vmin, vmax, nv = -1, 1, 7
Computes an interpolant of a 4D function via Hermite-tensor method.
ϵ = 0.1
γ = 3.0
Arguments:
- function_name : the name of the tested function.
Options_available: "f_3"
- nodes_types : an array of three strings with the spacing of nodes in corresponding directions.
Options available: :Uniform, :Chebyshev
- epsilon : a vector of the values of shape parameters in corresponding directions.
The values cannot be larger than 1!
- N : a vector with the numbers of basis functions in corresponding directions.
- Ne : a vector with the numbers of evaluation points in corresponding directions.
- interpolation_type : a string with the name of preferred basis.
Options available: :Hermite, :Radial
gamma - a factor allowing to scale the domain and the function accordingly.
hermite_x = Hermite(Chebyshev( xmin, xmax, nx ), ϵ, γ)
hermite_y = Hermite(Chebyshev( ymin, ymax, ny ), ϵ, γ)
hermite_z = Hermite(Chebyshev( zmin, zmax, nz ), ϵ, γ)
hermite_v = Hermite(Chebyshev( vmin, vmax, nv ), ϵ, γ)
© Anna Yurova, 2017
"""
function interpolate_4D(function_name, nodes_types, epsilon, N, Ne, interpolation_type, gamma)
xk = hermite_x.nodes
yk = hermite_y.nodes
zk = hermite_z.nodes
vk = hermite_v.nodes
# Set the nodes types in all directions.
nodes_type_x = nodes_types[1]
nodes_type_y = nodes_types[2]
nodes_type_z = nodes_types[3]
nodes_type_w = nodes_types[4]
nxe, nye, nze, nve = 8, 16, 16, 8
xe = collect(range(xmin, stop=xmax, length=nxe))
ye = collect(range(ymin, stop=ymax, length=nye))
ze = collect(range(zmin, stop=zmax, length=nze))
ve = collect(range(vmin, stop=vmax, length=nve))
dx = abs(xe[2] - xe[1])
dy = abs(ye[2] - ye[1])
dz = abs(ze[2] - ze[1])
dv = abs(ve[2] - ve[1])
## ----------------------- TEST CASES -------------------------------------
if function_name == :f_3
fk = [f(x, y, z, v) for x in xk, y in yk, z in zk, v in vk]
fe = [f(x, y, z, v) for x in xe, y in ye, z in ze, v in ve]
# Setup the interpolated function.
f = (x, y, z, w) -> cos(x.^2 + y.^2 + z.^2 + w.^2)
# Setup the bounds of the interpolation domain.
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
wmin = -1
wmax = 1
else
@error ("Wrong name of the test case! Possible options: f_3")
end
s = interpolate(hermite_x, hermite_y, hermite_z, hermite_v, fk, xe, ye, ze, ve)
## -------------- COLLOCATION AND EVALUATION MATRICES ---------------------
# Setup the basis used for the interpolation.
# Declare the functions for the generating of the collocation and evaluation
# matrices according to the chosed basis.
if interpolation_type == "Hermite"
generate_collocation_matrix_x = (xk, N, epsilon) -> evaluate_hermite(xk, N, epsilon, gamma)
generate_evaluation_matrix_x = (xk, xe, N, epsilon) -> evaluate_hermite(xe, N, epsilon, gamma)
max_error = maximum(abs.(s .- fe))
generate_collocation_matrix_y = (xk, N, epsilon) -> evaluate_hermite(xk, N, epsilon, gamma)
generate_evaluation_matrix_y = (xk, xe, N, epsilon) -> evaluate_hermite(xe, N, epsilon, gamma)
l2_error = sqrt(trapz((s - fe).^2, dx, dy, dz, dv))
generate_collocation_matrix_z = (zk, N, epsilon) -> evaluate_hermite(zk, N, epsilon, gamma)
generate_evaluation_matrix_z = (zk, ze, N, epsilon) -> evaluate_hermite(ze, N, epsilon, gamma)
println(" Hermite 4d max error : $max_error ")
println(" Hermite 4d L2 error : $l2_error ")
generate_collocation_matrix_w = (wk, N, epsilon) -> evaluate_hermite(wk, N, epsilon, gamma)
generate_evaluation_matrix_w = (wk, we, N, epsilon) -> evaluate_hermite(we, N, epsilon, gamma)
elseif interpolation_type == "Radial"
generate_collocation_matrix_x = (xk, N, epsilon) -> evaluate_radial(xk, xk, N, epsilon)
generate_evaluation_matrix_x = (xk, xe, N, epsilon) -> evaluate_radial(xk, xe, N, epsilon)
@test true
generate_collocation_matrix_y = (xk, N, epsilon) -> evaluate_radial(xk, xk, N, epsilon)
generate_evaluation_matrix_y = (xk, xe, N, epsilon) -> evaluate_radial(xk, xe, N, epsilon)
generate_collocation_matrix_z = (zk, N, epsilon) -> evaluate_radial(zk, zk, N, epsilon)
generate_evaluation_matrix_z = (zk, ze, N, epsilon) -> evaluate_radial(zk, ze, N, epsilon)
generate_collocation_matrix_w = (wk, N, epsilon) -> evaluate_radial(wk, wk, N, epsilon)
generate_evaluation_matrix_w = (wk, we, N, epsilon) -> evaluate_radial(wk, we, N, epsilon)
else
error("Wrong basis type! Possible options are: Radial, Hermite.")
end
## -------------------- NODES TYPES ---------------------------------------
# Declare the nodes generating functions in all directions according to the specified by the user.
# X - direction
if nodes_type_x == "Chebyshev"
generate_nodes_x = (Nx) -> generate_chebyshev_nodes(Nx, xmin, xmax)
elseif nodes_type_x == "Uniform"
generate_nodes_x = (Nx) -> generate_uniform_nodes(Nx, xmin, xmax)
else
error("Wrong node type! Possible options: Chebyshev, Uniform")