Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Allow arbitrary rectangular domain #47

Merged
merged 5 commits into from
Jun 13, 2023
Merged

Allow arbitrary rectangular domain #47

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
da72af1
rectangular code runs for square domain
JanJereczek Jun 13, 2023
56c2165
rect code runs 4 rect domain: u similar to square
JanJereczek Jun 13, 2023
a60b163
removed XMatrix
JanJereczek Jun 13, 2023
655cd3b
adapt test5 to rect domain
JanJereczek Jun 13, 2023
ada1e74
matrify() extended to examples
JanJereczek Jun 13, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 1 addition & 2 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ authors = ["Jan Swierczek-Jereczek"]
version = "0.1.0"

[deps]
AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Expand All @@ -16,7 +17,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"

[compat]
CUDA = "^4.0"
Expand All @@ -30,7 +30,6 @@ JLD2 = "^0.4"
OrdinaryDiffEq = "^6"
Reexport = "^1.2"
Statistics = "^1.9"
StatsBase = "~0.33"
julia = "^1.9"

[extras]
Expand Down
2 changes: 1 addition & 1 deletion docs/build/APIref.html
View file
Open in desktop

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion docs/build/search_index.js
View file
Open in desktop

Large diffs are not rendered by default.

3 changes: 1 addition & 2 deletions docs/src/APIref.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,8 +52,7 @@ update_slc_den!
years2seconds
seconds2years
m_per_sec2mm_per_yr
matrify_vectorconstant
matrify_constant
matrify
get_r
meshgrid
dist2angulardist
Expand Down
37 changes: 15 additions & 22 deletions src/FastIsostasy.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,25 +1,25 @@
module FastIsostasy

using LinearAlgebra
using Statistics: mean
using Distributions: MvNormal
using Interpolations: linear_interpolation, Flat
using FFTW
using StatsBase
using FastGaussQuadrature
using FFTW: fft, plan_fft, plan_ifft
using AbstractFFTs: Plan, ScaledPlan
using FastGaussQuadrature: gausslegendre
using DSP: conv
using CUDA
using CUDA: CuArray
using CUDA: CuArray, CUFFT.plan_fft, CUFFT.plan_ifft, allowscalar
using OrdinaryDiffEq: ODEProblem, solve, OrdinaryDiffEqAlgorithm
using EnsembleKalmanProcesses
using EnsembleKalmanProcesses.Observations
using EnsembleKalmanProcesses.ParameterDistributions

using Reexport
@reexport using Interpolations
@reexport using OrdinaryDiffEq: Euler, BS3, Tsit5, TanYam7, Vern9, VCABM, Rosenbrock23, QNDF, FBDF, ImplicitEuler
@reexport using OrdinaryDiffEq: Euler, Midpoint, Heun, Ralston, BS3, BS5, RK4, OwrenZen3,
OwrenZen4, OwrenZen5, Tsit5, DP5, RKO65, TanYam7, DP8, Feagin10, Feagin12, Feagin14,
TsitPap8, Vern6, Vern7, Vern8, Vern9, VCABM, Rosenbrock23, QNDF, FBDF, ImplicitEuler

# Euler, Midpoint, Heun, Ralston, RK4, OwrenZen3, OwrenZen4, OwrenZen5, DP5, RKO65,
# TanYam7, DP8, Feagin10, Feagin12, Feagin14, TsitPap8, BS5, Vern6, Vern7, Vern8,
# KuttaPRK2p5(dt=), Trapezoid(autodiff = false), PDIRK44(autodiff = false)

include("structs.jl")
Expand All @@ -30,7 +30,6 @@ include("mechanics.jl")
include("inversion.jl")

# structs.jl
export XMatrix
export ComputationDomain
export PhysicalConstants
export MultilayerEarth
Expand All @@ -39,30 +38,24 @@ export GeoState
export SuperStruct

# utils.jl
export years2seconds
export seconds2years
export m_per_sec2mm_per_yr
export years2seconds, seconds2years, m_per_sec2mm_per_yr
export meshgrid, dist2angulardist, latlon2stereo, stereo2latlon
export kernelpromote, convert2Array, copystructs2cpu

export meshgrid
export dist2angulardist
export latlon2stereo
export stereo2latlon

export kernelpromote
export convert2Array
export copystructs2cpu

export matrify_vectorconstant
export matrify, matrify
export loginterp_viscosity
export get_rigidity
export get_r
export gauss_distr

export samesize_conv

# derivatives.jl
export mixed_fdx
export mixed_fdy
export mixed_fdxx
export mixed_fdyy
export get_differential_fourier

# mechanics.jl
export quadrature1D
Expand Down
73 changes: 29 additions & 44 deletions src/derivatives.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,81 +1,81 @@

# FDM in x, 1st order
function central_fdx(M::XMatrix, h::T) where {T<:AbstractFloat}
function central_fdx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
n2 = size(M, 2)
return (view(M, :, 3:n2) - view(M, :, 1:n2-2)) ./ (2*h)
end

function forward_fdx(M::XMatrix, h::T) where {T<:AbstractFloat}
function forward_fdx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
return (view(M, :, 2) - view(M, :, 1)) ./ h
end

function backward_fdx(M::XMatrix, h::T) where {T<:AbstractFloat}
function backward_fdx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
n2 = size(M, 2)
return (view(M, :, n2) - view(M, :, n2-1)) ./ h
end

# FDM in y, 1st order
function mixed_fdx(M::XMatrix, h::T) where {T<:AbstractFloat}
function mixed_fdx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
return cat( forward_fdx(M,h), central_fdx(M,h), backward_fdx(M,h), dims=2 )
end

function central_fdy(M::XMatrix, h::T) where {T<:AbstractFloat}
function central_fdy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
n1 = size(M, 1)
return (view(M, 3:n1, :) - view(M, 1:n1-2, :)) ./ (2*h)
end

function forward_fdy(M::XMatrix, h::T) where {T<:AbstractFloat}
function forward_fdy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
return (view(M, 2, :) - view(M, 1, :)) ./ h
end

function backward_fdy(M::XMatrix, h::T) where {T<:AbstractFloat}
function backward_fdy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
n1 = size(M, 1)
return (view(M, n1, :) - view(M, n1-1, :)) ./ h
end

function mixed_fdy(M::XMatrix, h::T) where {T<:AbstractFloat}
function mixed_fdy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
return cat( forward_fdy(M,h)', central_fdy(M,h), backward_fdy(M,h)', dims=1 )
end

# FDM in x, 2nd order
function central_fdxx(M::XMatrix, h::T) where {T<:AbstractFloat}
function central_fdxx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
n2 = size(M, 2)
return (view(M, :, 3:n2) - 2 .* view(M, :, 2:n2-1) + view(M, :, 1:n2-2)) ./ h^2
end

function forward_fdxx(M::XMatrix, h::T) where {T<:AbstractFloat}
function forward_fdxx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
return (view(M, :, 3) - 2 .* view(M, :, 2) + view(M, :, 1)) ./ h^2
end

function backward_fdxx(M::XMatrix, h::T) where {T<:AbstractFloat}
function backward_fdxx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
n2 = size(M, 2)
return (view(M, :, n2) - 2 .* view(M, :, n2-1) + view(M, :, n2-2)) ./ h^2
end

function mixed_fdxx(M::XMatrix, h::T) where {T<:AbstractFloat}
function mixed_fdxx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
return cat( forward_fdxx(M,h), central_fdxx(M,h), backward_fdxx(M,h), dims=2 )
end

# FDM in y, 2nd order
function central_fdyy(M::XMatrix, h::T) where {T<:AbstractFloat}
function central_fdyy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
n1 = size(M, 1)
return (view(M, 3:n1, :) - 2 .* view(M, 2:n1-1, :) + view(M, 1:n1-2, :)) ./ h^2
end

function forward_fdyy(M::XMatrix, h::T) where {T<:AbstractFloat}
function forward_fdyy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
return (view(M, 3, :) - 2 .* view(M, 2, :) + view(M, 1, :)) ./ h^2
end

function backward_fdyy(M::XMatrix, h::T) where {T<:AbstractFloat}
function backward_fdyy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
n1 = size(M, 1)
return (view(M, n1, :) - 2 .* view(M, n1-1, :) + view(M, n1-2, :)) ./ h^2
end

function mixed_fdyy(M::XMatrix, h::T) where {T<:AbstractFloat}
function mixed_fdyy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat}
return cat( forward_fdyy(M,h)', central_fdyy(M,h), backward_fdyy(M,h)', dims=1 )
end

function mixed_fdxy(M::XMatrix, hx::T, hy::T) where {T<:AbstractFloat}
function mixed_fdxy(M::AbstractMatrix{T}, hx::T, hy::T) where {T<:AbstractFloat}
return mixed_fdy(mixed_fdx(M, hx), hy)
end

Expand All @@ -91,38 +91,23 @@ end

Compute the matrices representing the differential operators in the fourier space.
"""
function get_differential_fourier(
W::T,
N2::Int,
) where {T<:Real}
mu = T(π / W)
raw_coeffs = mu .* T.( vcat(0:N2, N2-1:-1:1) )
x_coeffs, y_coeffs = raw_coeffs, raw_coeffs
function get_differential_fourier(Wx::T, Wy::T, Nx::Int, Ny::Int) where {T<:Real}
mu_x = π / Wx
mu_y = π / Wy
x_coeffs = mu_x .* fftint(Nx)
y_coeffs = mu_y .* fftint(Ny)
X_coeffs, Y_coeffs = meshgrid(x_coeffs, y_coeffs)
harmonic_coeffs = X_coeffs .^ 2 + Y_coeffs .^ 2
pseudodiff_coeffs = sqrt.(harmonic_coeffs)
biharmonic_coeffs = harmonic_coeffs .^ 2
return pseudodiff_coeffs, harmonic_coeffs, biharmonic_coeffs
end

function precomp_fourier_dxdy(
M::XMatrix,
L1::T,
L2::T,
) where {T<:AbstractFloat}
n1, n2 = size(M)
k1 = 2 * π / L1 .* vcat(0:n1/2-1, 0, -n1/2+1:-1)
k2 = 2 * π / L2 .* vcat(0:n2/2-1, 0, -n2/2+1:-1)
p1, p2 = plan_fft(k1), plan_fft(k2)
ip1, ip2 = plan_ifft(k1), plan_ifft(k2)
return k1, k2, p1, p2, ip1, ip2
end

function fourier_dnx(
M::XMatrix,
k1::Vector{T},
p1::AbstractFFTs.Plan,
ip1::AbstractFFTs.ScaledPlan,
) where {T<:AbstractFloat}
return vcat( [real.(ip1 * ( ( im .* k1 ) .^ n .* (p1 * M[i, :]) ))' for i in axes(M,1)]... )
function fftint(N::Int)
N2 = N ÷ 2
if iseven(N)
return vcat(0:N2, N2-1:-1:1)
else
return vcat(0:N2, N2:-1:1)
end
end
13 changes: 5 additions & 8 deletions src/geostate.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,8 @@
Update the geoid of a `::GeoState` by convoluting the Green's function with the load anom.
"""
function update_geoid!(sstruct::SuperStruct{<:AbstractFloat})
sstruct.geostate.geoid .= view(
conv( sstruct.tools.geoidgreen, loadanom_green(sstruct) ),
sstruct.Omega.N2:2*sstruct.Omega.N-1-sstruct.Omega.N2,
sstruct.Omega.N2:2*sstruct.Omega.N-1-sstruct.Omega.N2,
)
sstruct.geostate.geoid .= samesize_conv(sstruct.tools.geoidgreen,
loadanom_green(sstruct), sstruct.Omega)
return nothing
end

Expand Down Expand Up @@ -109,12 +106,12 @@ end

"""

update_loadcolumns!(sstruct::SuperStruct, u::XMatrix, H_ice::XMatrix)
update_loadcolumns!(sstruct::SuperStruct, u::AbstractMatrix{T}, H_ice::AbstractMatrix{T})

Update the load columns of a `::GeoState`.
"""
function update_loadcolumns!(sstruct::SuperStruct{<:AbstractFloat},
u::XMatrix, H_ice::XMatrix)
function update_loadcolumns!(sstruct::SuperStruct{T}, u::AbstractMatrix{T},
H_ice::AbstractMatrix{T}) where {T<:AbstractFloat}

sstruct.geostate.b .= sstruct.refgeostate.b .+ u
sstruct.geostate.H_ice .= H_ice
Expand Down
Loading