From ed72c72929308c164421cc20b04008116d1b7772 Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Fri, 6 May 2022 16:27:33 +0200 Subject: [PATCH 01/10] initial commit kajiura non fft --- scripts/kajiura.jl | 51 ++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 47 insertions(+), 4 deletions(-) diff --git a/scripts/kajiura.jl b/scripts/kajiura.jl index 2f63720..ffba15b 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,42 @@ 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> 10.0 ) && (abs(d[y, x]) > 1e-8) + filter_nx_half = ceil(Int, n_h * h_max / Δx / 2) + filter_ny_half = 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 + 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,6 +136,17 @@ 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; 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(h_max, Δx, Δy, nx, ny, n_h) println(" Computing displacement matrix...") Threads.@threads for x ∈ 1:size(η, 2) @@ -191,6 +233,7 @@ function main() 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) + #apply_kajiura_fft!(b .+ current_disp, current_d_diff, current_η_diff, -maximum(b), -minimum(b), Δx, Δy) current_disp = d[:,:,t] println(" Writing output for timestep") From 60eee3f14241839fc20c63e35398e037fef6e369 Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Fri, 6 May 2022 17:05:04 +0200 Subject: [PATCH 02/10] fix nearcoast artifacts --- scripts/kajiura.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/scripts/kajiura.jl b/scripts/kajiura.jl index ffba15b..d5c4ce8 100644 --- a/scripts/kajiura.jl +++ b/scripts/kajiura.jl @@ -97,9 +97,10 @@ function apply_kajiura!(b::AbstractArray{Float64, 2}, d::AbstractArray{Float64, 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> 10.0 ) && (abs(d[y, x]) > 1e-8) - filter_nx_half = ceil(Int, n_h * h_max / Δx / 2) - filter_ny_half = ceil(Int, n_h * h_max / Δy / 2) + 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) From 9fb04bc18906f04a7200481c0e1a24b751449fcb Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Fri, 6 May 2022 17:29:24 +0200 Subject: [PATCH 03/10] fix fft implementation --- scripts/kajiura.jl | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/scripts/kajiura.jl b/scripts/kajiura.jl index d5c4ce8..7fbbd86 100644 --- a/scripts/kajiura.jl +++ b/scripts/kajiura.jl @@ -141,22 +141,21 @@ function compute_filter(h_max :: Float64, Δx :: Float64, Δy :: Float64, nx :: 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; water_level :: Float64 = 0., n_h :: Float64 = 20., - σ = precalculate_σ(h_min, h_max, Δx, Δy; n_h=n_h)) + 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(h_max, Δx, Δy, nx, ny, n_h) + 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) end @@ -178,8 +177,8 @@ function apply_kajiura_fft!(b::AbstractArray{Float64, 2}, d::AbstractArray{Float end function main() - if length(ARGS) != 3 - println("Usage: julia ./kajiura.jl ") + if length(ARGS) != 4 + println("Usage: julia ./kajiura.jl ") return end @@ -188,6 +187,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"] @@ -233,8 +240,12 @@ 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) - #apply_kajiura_fft!(b .+ current_disp, current_d_diff, current_η_diff, -maximum(b), -minimum(b), Δx, Δy) + + 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] println(" Writing output for timestep") From 5dcb8022c1ee2042dec78f74f443092722834aad Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Fri, 6 May 2022 17:44:21 +0200 Subject: [PATCH 04/10] update README --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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. From f07bab2db0ee6785fa216c432488d6767e2786af Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Wed, 18 May 2022 11:07:17 +0200 Subject: [PATCH 05/10] if no water, eta=d --- scripts/kajiura.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/scripts/kajiura.jl b/scripts/kajiura.jl index 7fbbd86..2d05490 100644 --- a/scripts/kajiura.jl +++ b/scripts/kajiura.jl @@ -115,6 +115,8 @@ function apply_kajiura!(b::AbstractArray{Float64, 2}, d::AbstractArray{Float64, 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 @@ -157,7 +159,7 @@ function apply_kajiura_fft!(b::AbstractArray{Float64, 2}, d::AbstractArray{Float if h_yx != 0. σ(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 From 2cf4a164e13b185d1a1ae5dcf280bed8a30e9ee5 Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Thu, 14 Jul 2022 10:44:44 +0200 Subject: [PATCH 06/10] locationFlag==2 should not be discarded --- src/rasterizer/rasterization.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/rasterizer/rasterization.jl b/src/rasterizer/rasterization.jl index 6512699..8dae83a 100644 --- a/src/rasterizer/rasterization.jl +++ b/src/rasterizer/rasterization.jl @@ -341,8 +341,9 @@ module Rasterization is_bath = is_bathy(m3n_simp_points, ctx.water_height) # locationFlag==1 is the acoustic side of an elastic-acoustic interface + # locationFlag==2 is an ordinary free surface boundary condition # https://seissol.readthedocs.io/en/latest/free-surface-output.html?#variables - if abs(ctx.locationFlag[simp_id]-1.0) < 1e-4 + if abs(ctx.locationFlag[simp_id]-1.0) < 1e-4 || abs(ctx.locationFlag[simp_id]-2.0) < 1e-4 l_bin_ids[simp_id] = ctx.n_threads * 2 + 1 # Unused bin, will be ignored later n_excluded_z_range += 1 continue From 00315014a55f32a452c652dd3a33f87d649558ee Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Mon, 18 Jul 2022 10:48:49 +0200 Subject: [PATCH 07/10] Revert "locationFlag==2 should not be discarded" This reverts commit 2cf4a164e13b185d1a1ae5dcf280bed8a30e9ee5. --- src/rasterizer/rasterization.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/rasterizer/rasterization.jl b/src/rasterizer/rasterization.jl index 8dae83a..6512699 100644 --- a/src/rasterizer/rasterization.jl +++ b/src/rasterizer/rasterization.jl @@ -341,9 +341,8 @@ module Rasterization is_bath = is_bathy(m3n_simp_points, ctx.water_height) # locationFlag==1 is the acoustic side of an elastic-acoustic interface - # locationFlag==2 is an ordinary free surface boundary condition # https://seissol.readthedocs.io/en/latest/free-surface-output.html?#variables - if abs(ctx.locationFlag[simp_id]-1.0) < 1e-4 || abs(ctx.locationFlag[simp_id]-2.0) < 1e-4 + if abs(ctx.locationFlag[simp_id]-1.0) < 1e-4 l_bin_ids[simp_id] = ctx.n_threads * 2 + 1 # Unused bin, will be ignored later n_excluded_z_range += 1 continue From ebcf7111efa5d3e0ba56faf13fa215238d81283c Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Fri, 14 Jun 2024 08:48:25 +0200 Subject: [PATCH 08/10] fix when no locationFlag --- src/rasterizer/rasterization.jl | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) 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 ,)) From 5aa1f7b5474674625f275b2c043b4df76d0c941b Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Fri, 14 Jun 2024 10:28:07 +0200 Subject: [PATCH 09/10] change outputed variables in kajiura filter --- scripts/kajiura.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/scripts/kajiura.jl b/scripts/kajiura.jl index 2d05490..243c3b7 100644 --- a/scripts/kajiura.jl +++ b/scripts/kajiura.jl @@ -222,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)) @@ -232,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") @@ -249,10 +250,10 @@ function main() 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 From 2848cc00f9ce3cd74d0cbfdf74f5176104139053 Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Thu, 25 Jul 2024 09:37:45 +0200 Subject: [PATCH 10/10] fix #13 (order of netcdf dims) --- src/io/netcdf.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) 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