Commit 04360d71 authored by Pierre Navaro's avatar Pierre Navaro
Browse files

add nodes and interpolation types in src

parent 05ca882b
This diff is collapsed.
module HermiteGF
using Distributed
using SharedArrays
export interpolate_1D
export interpolate_2D
export Hermite, Radial
export Chebyshev, Uniform
export interpolate
include("helper_functions.jl")
include("ndgrid.jl")
include("evaluate_tensor_product.jl")
#@everywhere include("evaluate_tensor_product_parallel.jl")
include("evaluate_hermite.jl")
" Interpolation type (Hermite or Radial) "
abstract type InterpolationType end
" Node positions (Uniform or Chebyshev) "
abstract type NodesType end
include("chebyshev.jl")
include("uniform.jl")
include("hermite.jl")
include("radial.jl")
include("interpolation1d.jl")
include("interpolation2d.jl")
include("interpolation3d.jl")
......
"""
Chebyshev( nx, xmin, xmax )
Chebyshev nodes
"""
struct Chebyshev <: NodesType
nx :: Int64
xmin :: Float64
xmax :: Float64
xk :: Vector{Float64}
function Chebyshev( xmin, xmax, nx )
xk = zeros(Float64, nx)
θ = (π:-π/(nx-1):-1e-14);
xk .= ((cos.(θ).+1)*(xmax-xmin))/2 .+ xmin
new( nx, xmin, xmax, xk )
end
end
using TensorOperations
"""
evaluate_s(X_all, Y_all, Z_all, W_all, F)
evaluate_s(x, y, z, v, w, f)
In this file the functions for the serial computing of the interpolated values 's'
in 1-4D are implemented. Here the tensor representation of 's' is used (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, v, w, f)
Computing s = (W_all ⊗ Z_all ⊗ Y_all ⊗ X_all) vec(F) (see Sec. 4.1 of the paper)
# 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)
@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]
end
s
end
"""
evaluate_s(x, y, z, w, f)
Computing s = (w ⊗ z ⊗ y ⊗ x) vec(f) (see Sec. 4.1 of the paper)
"""
function evaluate_s(X_all, Y_all, Z_all, W_all, F)
function evaluate_s(x, y, z, w, f)
# Extract the number of collocation and evaluation points in each dimension
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
Ny = size(Y_all, 2);
Ne2 = size(Y_all, 1);
Nz = size(Z_all, 2);
Ne3 = size(Z_all, 1);
Nw = size(W_all, 2);
Ne4 = size(W_all, 1);
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);
s = zeros(ne1, ne2, ne3, ne4)
@inbounds for col_dim_4 = 1:Nw
@inbounds for eval_dim_4 = 1:Ne4
@inbounds for col_dim_3 = 1:Nz
@inbounds for eval_dim_3 = 1:Ne3
@inbounds for col_dim_2 = 1:Ny
@inbounds for eval_dim_2 = 1:Ne2
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1, eval_dim_2, eval_dim_3, eval_dim_4] = (
s[eval_dim_1,eval_dim_2, eval_dim_3, eval_dim_4]
+ X_all[eval_dim_1, col_dim_1] * Y_all[eval_dim_2, col_dim_2]
* Z_all[eval_dim_3, col_dim_3] * W_all[eval_dim_4, col_dim_4]
* F[col_dim_1, col_dim_2, col_dim_3, col_dim_4])
end
end
end
end
end
end
end
@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]
end
s
end
"""
evaluate_s(X_all, Y_all, Z_all, F)
evaluate_s(x, y, z, f)
Extract the number of collocation and evaluation points in each dimension
Computing `s = (Z_all ⊗ Y_all ⊗ X_all) 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_all, Y_all, Z_all, F)
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
Ny = size(Y_all, 2);
Ne2 = size(Y_all, 1);
Nz = size(Z_all, 2);
Ne3 = size(Z_all, 1);
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);
s = zeros(ne1,ne2,ne3)
@inbounds for col_dim_3 = 1:Nz
@inbounds for eval_dim_3 = 1:Ne3
@inbounds for col_dim_2 = 1:Ny
@inbounds for eval_dim_2 = 1:Ne2
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1, eval_dim_2, eval_dim_3] = (
s[eval_dim_1,eval_dim_2, eval_dim_3]
+ X_all[eval_dim_1, col_dim_1]
* Y_all[eval_dim_2, col_dim_2]
* Z_all[eval_dim_3 ,col_dim_3]
* F[col_dim_1, col_dim_2, col_dim_3] )
end
end
end
end
end
@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
......@@ -89,55 +72,40 @@ end
"""
evaluate_s(X_all, Y_all, F)
evaluate_s(x, y, f)
Extract the number of collocation and evaluation points in each dimension
Computing `s = (Y_all ⊗ X_all) vec(F)` (see Sec. 4.1 of the paper)
Computing `s = (y ⊗ x) vec(f)` (see Sec. 4.1 of the paper)
"""
function evaluate_s(X_all, Y_all, F)
function evaluate_s(x, y, f)
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
Ny = size(Y_all, 2);
Ne2 = size(Y_all, 1);
ne1 = size(x)[1]
ne2 = size(y)[1]
# Initialize the result tensor
s = zeros(Ne1, Ne2);
@inbounds for col_dim_2 = 1:Ny
@inbounds for eval_dim_2 = 1:Ne2
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1, eval_dim_2] = ( s[eval_dim_1,eval_dim_2]
+ X_all[eval_dim_1, col_dim_1]
* Y_all[eval_dim_2, col_dim_2]*F[col_dim_1, col_dim_2])
end
end
end
s = zeros(ne1, ne2)
@tensor begin
s[e1,e2] = s[e1,e2]+x[e1,c1]*y[e2,c2]*f[c1,c2]
end
s
end
"""
evaluate_s(X_all, F)
evaluate_s(x, f)
Extract the number of collocation and evaluation points
Computing s = X_all * vec(F) (see Sec. 4.1 of the paper)
Computing s = x * vec(f) (see Sec. 4.1 of the paper)
"""
function evaluate_s(X_all, F)
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
function evaluate_s(x, f)
# Initialize the result tensor
s = zeros(Ne1);
s = similar(x)
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1] = s[eval_dim_1] + X_all[eval_dim_1, col_dim_1]*F[col_dim_1];
end
@tensor begin
s[e1] = s[e1] + x[e1,c1]*f[c1]
end
s
......
using SharedArrays
###
# In this file the functions for the parallel computing of the interpolated values 's'
# in 4 and 5D are implemented. Here the tensor representation of 's' is used (see Sec. 4.1 of the paper)
###
"""
This function is what we call from the rest of the code.
It distributes the work among processes and collects the results.
"""
function evaluate_s_4D_exec!(X_all, Y_all, Z_all, W_all, F, s)
@sync begin
for p in procs(s)
@show p
@async remotecall_wait(evaluate_s_4D_parallel_chunk!, p, X_all, Y_all, Z_all, W_all, F, s)
end
end
return s;
end
"""
This function is what we call from the rest of the code.
It distributes the work among processes and collects the results.
"""
function evaluate_s_5D_exec!(X_all, Y_all, Z_all, W_all, V_all, F, s)
@sync begin
for p in procs(s)
@show p
@async remotecall_wait(evaluate_s_5D_parallel_chunk!, p, X_all, Y_all, Z_all, W_all, V_all, F, s)
end
end
return s;
end
"""
This function retuns the (irange,jrange) indexes assigned to a worker.
In the case of the HermiteGF code, q corresponds to the tensor 's' containing
the values of the interpolant at the evaluation points
"""
function myrange_4D(q::SharedArray)
idx = indexpids(q)
if idx == 0
# This worker is not assigned a piece
return 1:0
end
# Split the number of points in 4-th dimension into <number of threads>
# equal chuncks
nchunks = length(procs(q))
splits = [round(Int, s) for s in linspace(0,size(q,4),nchunks+1)]
return splits[idx]+1:splits[idx+1]
end
"""
This function retuns the (irange,jrange) indexes assigned to a worker
In the case of the HermiteGF code, q corresponds to the tensor 's' containing
the values of the interpolant at the evaluation points
"""
function myrange_5D(q::SharedArray)
idx = indexpids(q)
if idx == 0
# This worker is not assigned a piece
return 1:0
end
# Split the number of points in 5-th dimension into <number of threads>
# equal chuncks
nchunks = length(procs(q))
splits = [round(Int, s) for s in linspace(0,size(q,5),nchunks+1)]
return splits[idx]+1:splits[idx+1]
end
"""
This is the kernel where the actual multiplication happens.
Every thread takes care of the slices within 'irange' of the evaluation points in 4-th dimension
"""
function evaluate_s_4D_parallel!(X_all, Y_all, Z_all, W_all, F, s, irange)
#@show (irange) # display what irange each worker is taking care of
# Extract the number of collocation and evaluation points in each dimension
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
Ny = size(Y_all, 2);
Ne2 = size(Y_all, 1);
Nz = size(Z_all, 2);
Ne3 = size(Z_all, 1);
Nw = size(W_all, 2);
Ne4 = size(W_all, 1);
# Computing s(:, :, :, irange) = (W_all(irange, :) \otimes Z_all \otimes Y_all \otimes X_all) vec(F) (see Sec. 4.1 of the paper)
for eval_dim_4 in irange
@inbounds for col_dim_4 = 1:Nw
@inbounds for col_dim_3 = 1:Nz
@inbounds for eval_dim_3 = 1:Ne3
@inbounds for col_dim_2 = 1:Ny
@inbounds for eval_dim_2 = 1:Ne2
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1, eval_dim_2, eval_dim_3, eval_dim_4] = (
s[eval_dim_1,eval_dim_2, eval_dim_3, eval_dim_4]
+ X_all[eval_dim_1, col_dim_1]
* Y_all[eval_dim_2, col_dim_2]
* Z_all[eval_dim_3, col_dim_3]
* W_all[eval_dim_4, col_dim_4]
* F[col_dim_1, col_dim_2, col_dim_3, col_dim_4])
end
end
end
end
end
end
end
end
s
end
"""
This is the kernel where the actual multiplication happens.
Every thread takes care of the slices within the 'irange' of the evaluation points in 5-th dimension
"""
function evaluate_s_5D_parallel!(X_all, Y_all, Z_all, W_all, V_all, F, s, irange)
#@show (irange) # display what irange each worker is taking care of
# Extract the number of collocation and evaluation points in each dimension
Nx = size(X_all, 2);
Ne1 = size(X_all, 1);
Ny = size(Y_all, 2);
Ne2 = size(Y_all, 1);
Nz = size(Z_all, 2);
Ne3 = size(Z_all, 1);
Nw = size(W_all, 2);
Ne4 = size(W_all, 1);
Nv = size(V_all, 2);
Ne5 = size(V_all, 1);
# Computing s(:, :, :, :, irange) = (V_all(irange, :) \otimes W_all \otimes Z_all \otimes Y_all \otimes X_all) vec(F) (see Sec. 4.1 of the paper)
for eval_dim_5 in irange
@inbounds for col_dim_5 = 1:Nv
@inbounds for col_dim_4 = 1:Nw
@inbounds for eval_dim_4 = 1:Ne4
@inbounds for col_dim_3 = 1:Nz
@inbounds for eval_dim_3 = 1:Ne3
@inbounds for col_dim_2 = 1:Ny
@inbounds for eval_dim_2 = 1:Ne2
@inbounds for col_dim_1 = 1:Nx
@inbounds for eval_dim_1 = 1:Ne1
s[eval_dim_1, eval_dim_2, eval_dim_3, eval_dim_4, eval_dim_5] = (
s[eval_dim_1,eval_dim_2, eval_dim_3, eval_dim_4, eval_dim_5]
+ X_all[eval_dim_1, col_dim_1]
* Y_all[eval_dim_2, col_dim_2]
* Z_all[eval_dim_3, col_dim_3]
* W_all[eval_dim_4, col_dim_4]
* V_all[eval_dim_5, col_dim_5]
* F[col_dim_1, col_dim_2, col_dim_3, col_dim_4, col_dim_5])
end
end
end
end
end
end
end
end
end
end
s
end
# This is an interface simplifying the call of the function from the rest of the code
evaluate_s_4D_parallel_chunk!(X_all, Y_all, Z_all, W_all, F, s) = evaluate_s_4D_parallel!(X_all, Y_all, Z_all, W_all, F, s, myrange_4D(s))
evaluate_s_5D_parallel_chunk!(X_all, Y_all, Z_all, W_all, V_all, F, s) = evaluate_s_5D_parallel!(X_all, Y_all, Z_all, W_all, V_all, F, s, myrange_5D(s))
"""
Hermite( nodes::NodesType, epsilon, gamma )
Hermite inteprolation
"""
mutable struct Hermite <: InterpolationType
nx :: Vector{Int64}
epsilon :: Vector{Float64}
gamma :: Vector{Float64}
xk :: Vector{Array{Float64,1}}
colloc_mat :: AbstractVecOrMat{Float64}
function Hermite( nodes::Vector{NodesType},
epsilon::Vector{Float64},
gamma::Vector{Float64} )
xk = nodes.xk
nx = nodes.nx
colloc_mat = evaluate_hermite(xk, nx, epsilon, gamma)
new( nx, epsilon, gamma, xk, colloc_mat )
end
end
"""
Hermite( xe )
returns evaluation matrix
"""
function (interp::Hermite)( xe::Array{Float64,1} )
nx = interp.nx
epsilon = interp.epsilon
gamma = interp.gamma
evaluate_hermite( xe, nx, epsilon, gamma)
end
function interpolate(interp::InterpolationType,
f::Array{Float64,1}, xe::Array{Float64,1} )
ϵ = interp.epsilon
nx = interp.nx
xk = interp.xk
# Precompute 1D interpolation matrix
x_all = interp(xe) / interp.colloc_mat
# Compute the interpolated values
ne = Int(length(xe))
# Initialize the result tensor
s = similar(xe)
@tensor begin
s[e] = s[e] + x_all[e, c]*f[c]
end
s
end
"""
Radial( nodes::NodesType, epsilon, gamma )
Radial inteprolation
"""
mutable struct Radial <: InterpolationType
nx :: Int64
epsilon :: Float64
xk :: Array{Float64,1}
colloc_mat :: Array{Float64,2}
function Radial( nodes::NodesType, epsilon )
xk = nodes.xk
nx = nodes.nx
colloc_mat = evaluate_radial(xk, xk, nx, epsilon)
new( nx, epsilon, xk, colloc_mat )
end
end
"""
Radial( xe )
returns evaluation matrix
"""
function (interp::Radial)( xe::Array{Float64,1} )
nx = interp.nx
epsilon = interp.epsilon
xk = interp.xk
evaluate_radial(xk, xe, nx, epsilon)
end
"""
Uniform( nx, xmin, xmax )
Uniform nodes
"""
struct Uniform <: NodesType
nx :: Int64
xmin :: Float64
xmax :: Float64
xk :: Vector{Float64}
function Uniform( xmin, xmax, nx )
xk = zeros(Float64, nx)
xk .= range(-1, stop=1, length=nx);
xk .= ((xk .+ 1) * (xmax - xmin)) /2 .+ xmin
new( nx, xmin, xmax, xk )
end
end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment