Thomas/kajiura filter nofft

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 <in_file.nc> <out_file.nc> <timestep_end>
+julia --project scripts/kajiura.jl <in_file.nc> <out_file.nc> <timestep_end> <filter_depth (0: variable, slower code)>
 ```
 
 _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.

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 <in_file.nc> <out_file.nc> <timestep_end>")
+    if length(ARGS) != 4
+        println("Usage: julia ./kajiura.jl <in_file.nc> <out_file.nc> <timestep_end> <filter_depth (0: variable, slower code)>")
         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
 

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
 
 

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 ,))