diff --git a/README.md b/README.md index 2748332..1b4899e 100644 --- a/README.md +++ b/README.md @@ -59,7 +59,7 @@ More help for running on the cluster can be found [here][5]. Once you have rasterized the SeisSol output files, you can optionally use Kajiura's filter on the outputs: ```bash -julia --project scripts/kajiura.jl +julia --project scripts/kajiura.jl ``` _Note: Kajiura's Filter needs to process all timesteps prior to `timestep_end`._ @@ -82,7 +82,7 @@ julia --project src/main.jl -m 8G --water-height=2000 -o ~/sampler-output.nc --s Only rasterizes the 2D _seafloor_ elevation over time. Optionally thereafter: ```bash -julia --project src/main.jl ~/sampler-output.nc ~/kajiura-output.nc 300 +julia --project src/main.jl ~/sampler-output.nc ~/kajiura-output.nc 300 0. ``` Applies Kajiura's Filter for 300 timesteps. diff --git a/scripts/kajiura.jl b/scripts/kajiura.jl index 2f63720..243c3b7 100644 --- a/scripts/kajiura.jl +++ b/scripts/kajiura.jl @@ -35,7 +35,7 @@ const G = begin end function precalculate_σ(h_min :: Float64, h_max :: Float64, Δx :: Float64, Δy :: Float64; n_h :: Float64 = 20.) - Δh = (h_max - h_min) / 10 + Δh = (h_max - max(0.0, h_min)) / 10 if Δh == 0.; Δh = 1.; end l_h = max(0.0001, h_min):Δh:(h_max + Δh) @@ -88,11 +88,45 @@ function apply_kajiura!(b::AbstractArray{Float64, 2}, d::AbstractArray{Float64, h_min :: Float64, h_max :: Float64, Δx :: Float64, Δy :: Float64; water_level :: Float64 = 0., n_h :: Float64 = 20., σ = precalculate_σ(h_min, h_max, Δx, Δy; n_h=n_h)) - filter_nx_half = ceil(Int, n_h * h_max / Δx / 2) - filter_ny_half = ceil(Int, n_h * h_max / Δy / 2) - nx = size(η, 2) ny = size(η, 1) + filter_nx_half0 = ceil(Int, n_h * h_max / Δx / 2) + filter_ny_half0 = ceil(Int, n_h * h_max / Δy / 2) + + println(" Computing displacement matrix...") + Threads.@threads for x ∈ 1:nx + for y ∈ 1:ny + h_yx = max(0., water_level - b[y, x]) # height 0 on land + if (h_yx> 0.01 ) && (abs(d[y, x]) > 1e-8) + # the min is required to avoid aritacts at very shallow depth + filter_nx_half = min(3, ceil(Int, n_h * h_max / Δx / 2)) + filter_ny_half = min(3, ceil(Int, n_h * h_max / Δy / 2)) + + xmin = max(1, x-filter_nx_half) + xmax = min(nx, x+filter_nx_half) + + ymin = max(1, y-filter_ny_half) + ymax = min(ny, y+filter_ny_half) + + filter = zeros(Float64, (ymax-ymin+1, xmax-xmin+1)) + for x1 ∈ xmin:xmax + for y1 ∈ ymin:ymax + filter[y1-ymin+1,x1-xmin+1] = G(sqrt(((x1-x) * Δx)^2 + ((y1-y) * Δy)^2) / h_yx) + end + end + η[ymin:ymax, xmin:xmax] += σ(h_yx) * Δx * Δy / h_yx^2 * d[y, x] * filter + else + η[y,x] = d[y, x] + end + end + end +end + +function compute_filter(h_max :: Float64, Δx :: Float64, Δy :: Float64, nx :: Int64, ny :: Int64, + n_h :: Float64) + + filter_nx_half = ceil(Int, n_h * h_max / Δx / 2) + filter_ny_half = ceil(Int, n_h * h_max / Δy / 2) filter = zeros(Float64, (ny, nx)) @@ -105,17 +139,27 @@ function apply_kajiura!(b::AbstractArray{Float64, 2}, d::AbstractArray{Float64, filter[((ny + y) % ny) + 1, ((nx + x) % nx) + 1] = G(sqrt((x * Δx)^2 + (y * Δy)^2) / h_max) end end + return filter +end + +function apply_kajiura_fft!(b::AbstractArray{Float64, 2}, d::AbstractArray{Float64, 2}, η::AbstractArray{Float64, 2}, + h_min :: Float64, h_max :: Float64, Δx :: Float64, Δy :: Float64, filter_depth :: Float64; + water_level :: Float64 = 0., n_h :: Float64 = 20., σ = precalculate_σ(h_min, h_max, Δx, Δy; n_h=n_h)) + + nx = size(η, 2) + ny = size(η, 1) + + filter = compute_filter(filter_depth, Δx, Δy, nx, ny, n_h) println(" Computing displacement matrix...") Threads.@threads for x ∈ 1:size(η, 2) for y ∈ 1:size(η, 1) h_yx = max(0., water_level - b[y, x]) # height 0 on land - η[y, x] = if h_yx != 0. - σ(h_yx) * Δx * Δy / h_yx^2 * d[y, x] + σ(filter_depth) * Δx * Δy / filter_depth^2 * d[y, x] else - 0. # No displacement where no water is (Kajiura not applicable) + d[y, x] # No displacement where no water is (Kajiura not applicable) end end end @@ -135,8 +179,8 @@ function apply_kajiura!(b::AbstractArray{Float64, 2}, d::AbstractArray{Float64, end function main() - if length(ARGS) != 3 - println("Usage: julia ./kajiura.jl ") + if length(ARGS) != 4 + println("Usage: julia ./kajiura.jl ") return end @@ -145,6 +189,14 @@ function main() in_filename = ARGS[1] out_filename = ARGS[2] t_end = parse(Int, ARGS[3]) + filter_depth = parse(Float64, ARGS[4]) + + if filter_depth!=0. + println("Using fft implementation with filter_depth: $(filter_depth).") + else + println("Using slow implementation with variable seafloor depth.") + end + nc = NetCDF.open(in_filename, mode=NC_NOWRITE) d = nc["d"] @@ -170,6 +222,7 @@ function main() b = Array{Float64, 2}(undef, (ny, nx)) ncread!(in_filename, "b", b) current_disp = zeros(ny, nx) + current_eta = zeros(ny, nx) current_η_diff = Array{Float64, 2}(undef, (ny, nx)) current_d_diff = Array{Float64, 2}(undef, (ny, nx)) @@ -180,8 +233,8 @@ function main() xatts = Dict("units" => "m") yatts = Dict("units" => "m") - nccreate(out_filename, "eta_diff", "y", ly, yatts, "x", lx, xatts, "time", l_times[1:t_end+1], timeatts) - nccreate(out_filename, "d_diff", "y", "x", "time") + nccreate(out_filename, "eta", "y", ly, yatts, "x", lx, xatts, "time", l_times[1:t_end+1], timeatts) + nccreate(out_filename, "d", "y", "x", "time") nccreate(out_filename, "b", "y", "x") ncwrite(b, out_filename, "b") @@ -190,12 +243,17 @@ function main() println("Working on timestep $(t - 1) of $(t_end)") current_η_diff .= 0. current_d_diff .= d[:,:,t] .- current_disp - apply_kajiura!(b .+ current_disp, current_d_diff, current_η_diff, -maximum(b), -minimum(b), Δx, Δy) - current_disp = d[:,:,t] + if filter_depth!=0. + apply_kajiura_fft!(b .+ current_disp, current_d_diff, current_η_diff, -maximum(b), -minimum(b), Δx, Δy, filter_depth) + else + apply_kajiura!(b .+ current_disp, current_d_diff, current_η_diff, -maximum(b), -minimum(b), Δx, Δy) + end + current_disp = d[:,:,t] + current_eta = current_eta[:,:] + current_η_diff println(" Writing output for timestep") - ncwrite(current_d_diff, out_filename, "d_diff", start=[1,1,t], count=[-1,-1,1]) - ncwrite(current_η_diff, out_filename, "eta_diff", start=[1,1,t], count=[-1,-1,1]) + ncwrite(current_disp, out_filename, "d", start=[1,1,t], count=[-1,-1,1]) + ncwrite(current_eta, out_filename, "eta", start=[1,1,t], count=[-1,-1,1]) end end diff --git a/src/io/netcdf.jl b/src/io/netcdf.jl index 3fa8791..95a741e 100644 --- a/src/io/netcdf.jl +++ b/src/io/netcdf.jl @@ -24,7 +24,7 @@ module NC For each key in those var mappings, a variable corresponding to the mapped output name of that key is created, and assigned the attributes (e.g. units) of the input variable. """ - function create_netcdf(filename :: AbstractString, x_vals, y_vals, t_vals, + function create_netcdf(filename :: AbstractString, y_vals, x_vals, t_vals, static_var_mappings :: Array{Args.VarMapping, 1}, dynamic_var_mappings :: Array{Args.VarMapping, 1}) :: NcFile @@ -51,13 +51,12 @@ module NC static_vars = [] dynamic_vars = [] - for stat_mapping ∈ static_var_mappings - append!(static_vars, map(in_var -> NcVar(stat_mapping[in_var], [y_dim, x_dim], atts=get_units(in_var)), collect(keys(stat_mapping)))) + append!(static_vars, map(in_var -> NcVar(stat_mapping[in_var], [x_dim, y_dim], atts=get_units(in_var)), collect(keys(stat_mapping)))) end for dyn_mapping ∈ dynamic_var_mappings - append!(dynamic_vars, map(in_var -> NcVar(dyn_mapping[in_var], [y_dim, x_dim, t_dim], atts=get_units(in_var)), collect(keys(dyn_mapping)))) + append!(dynamic_vars, map(in_var -> NcVar(dyn_mapping[in_var], [x_dim, y_dim, t_dim], atts=get_units(in_var)), collect(keys(dyn_mapping)))) end diff --git a/src/rasterizer/rasterization.jl b/src/rasterizer/rasterization.jl index 6512699..712dbc9 100644 --- a/src/rasterizer/rasterization.jl +++ b/src/rasterizer/rasterization.jl @@ -224,8 +224,19 @@ module Rasterization simplices, points = grid_of(xdmf) n_dims = size(simplices, 1) - 1 + n_elements = size(simplices, 2) if n_dims==2 - locationFlag = data_of(xdmf, 1 , "locationFlag") + try + locationFlag = data_of(xdmf, 1, "locationFlag") + catch e + if isa(e, ErrorException) + println("locationFlag not found, assuming locationFlag=2...") + locationFlag = fill(2.0, (n_elements,)) + else + println("Unhandled exception: ", e) + rethrow(e) # Rethrow the error if it's not a LoadError + end + end else #for volume output we use a dummy array adhering the locationFlag format in the structure RasterizationContext locationFlag = fill(-1.0, (1 ,))