`GaussianMask` fails for domains with `Flat` topology

[xkykai]

See the MWE:

using Oceananigans

grid = RectilinearGrid(CPU(), Float64,
                        size = (8, 8), 
                        x = (0, 1),
                        z = (0, 1),
                        topology = (Bounded, Flat, Bounded))

mask = GaussianMask{:z}(center=grid.Lz, width=grid.Lz/10)

uvw_forcing = Relaxation(rate=0.1, mask=mask)

model = NonhydrostaticModel(; grid, forcing = (u=uvw_forcing, v=uvw_forcing, w=uvw_forcing))

simulation = Simulation(model; Δt=0.1, stop_time=1)

run!(simulation)

which gives an error of

julia> run!(simulation)
[ Info: Initializing simulation...
ERROR: TaskFailedException

    nested task error: MethodError: no method matching (::GaussianMask{:z, Float64})(::Float64, ::Float64)
    The object of type `GaussianMask{:z, Float64}` exists, but no method is defined for this combination of argument types when trying to treat it as a callable object.

    Closest candidates are:
      (::GaussianMask{:z})(::Any, ::Any, ::Any)
       @ Oceananigans c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Forcings\relaxation.jl:151

    Stacktrace:
     [1] Relaxation
       @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Forcings\relaxation.jl:90 [inlined]
     [2] ContinuousForcing
       @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Forcings\continuous_forcing.jl:141 [inlined]
     [3] u_velocity_tendency
       @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Models\NonhydrostaticModels\nonhydrostatic_tendency_kernel_functions.jl:93 [inlined]
     [4] cpu_compute_Gu!
       @ C:\Users\xinle\.julia\packages\KernelAbstractions\X5fk1\src\macros.jl:314 [inlined]
     [5] cpu_compute_Gu!(__ctx__::KernelAbstractions.CompilerMetadata{…}, Gu::Field{…}, grid::RectilinearGrid{…}, args::Tuple{…})
       @ Oceananigans.Models.NonhydrostaticModels .\none:0
     [6] __thread_run(tid::Int64, len::Int64, rem::Int64, obj::KernelAbstractions.Kernel{…}, ndrange::Nothing, iterspace::KernelAbstractions.NDIteration.NDRange{…}, args::Tuple{…}, dynamic::KernelAbstractions.NDIteration.NoDynamicCheck)
       @ KernelAbstractions C:\Users\xinle\.julia\packages\KernelAbstractions\X5fk1\src\cpu.jl:145
     [7] (::KernelAbstractions.var"#29#32"{…})()
       @ KernelAbstractions C:\Users\xinle\.julia\packages\KernelAbstractions\X5fk1\src\cpu.jl:120

...and 1 more exception.

Stacktrace:
  [1] sync_end(c::Channel{Any})
    @ Base .\task.jl:466
  [2] macro expansion
    @ .\task.jl:499 [inlined]
  [3] __run(obj::KernelAbstractions.Kernel{…}, ndrange::Nothing, iterspace::KernelAbstractions.NDIteration.NDRange{…}, args::Tuple{…}, dynamic::KernelAbstractions.NDIteration.NoDynamicCheck, static_threads::Bool)
    @ KernelAbstractions C:\Users\xinle\.julia\packages\KernelAbstractions\X5fk1\src\cpu.jl:119
  [4] (::KernelAbstractions.Kernel{…})(::Field{…}, ::Vararg{…}; ndrange::Nothing, workgroupsize::Nothing)
    @ KernelAbstractions C:\Users\xinle\.julia\packages\KernelAbstractions\X5fk1\src\cpu.jl:46
  [5] (::KernelAbstractions.Kernel{…})(::Field{…}, ::Vararg{…})
    @ KernelAbstractions C:\Users\xinle\.julia\packages\KernelAbstractions\X5fk1\src\cpu.jl:39
  [6] _launch!(::CPU, ::RectilinearGrid{…}, ::Symbol, ::Function, ::Field{…}, ::RectilinearGrid{…}, ::Vararg{…}; exclude_periphery::Bool, reduced_dimensions::Tuple{}, active_cells_map::Nothing)
    @ Oceananigans.Utils c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Utils\kernel_launching.jl:378
  [7] _launch!
    @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Utils\kernel_launching.jl:365 [inlined]
  [8] launch!
    @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Utils\kernel_launching.jl:342 [inlined]
  [9] #compute_interior_tendency_contributions!#25
    @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Models\NonhydrostaticModels\compute_nonhydrostatic_tendencies.jl:97 [inlined]
 [10] compute_interior_tendency_contributions!
    @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Models\NonhydrostaticModels\compute_nonhydrostatic_tendencies.jl:49 [inlined]
 [11] compute_tendencies!(model::NonhydrostaticModel{…}, callbacks::Vector{…})
    @ Oceananigans.Models.NonhydrostaticModels c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Models\NonhydrostaticModels\compute_nonhydrostatic_tendencies.jl:36
 [12] #apply_regionally!#62
    @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Utils\multi_region_transformation.jl:88 [inlined]
 [13] apply_regionally!
    @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Utils\multi_region_transformation.jl:85 [inlined]
 [14] macro expansion
    @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Utils\multi_region_transformation.jl:146 [inlined]
 [15] update_state!(model::NonhydrostaticModel{…}, callbacks::Vector{…}; compute_tendencies::Bool)
    @ Oceananigans.Models.NonhydrostaticModels c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Models\NonhydrostaticModels\update_nonhydrostatic_model_state.jl:53
 [16] update_state! (repeats 2 times)
    @ c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Models\NonhydrostaticModels\update_nonhydrostatic_model_state.jl:21 [inlined]
 [17] initialize!(sim::Simulation{…})
    @ Oceananigans.Simulations c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Simulations\run.jl:210
 [18] time_step!(sim::Simulation{…})
    @ Oceananigans.Simulations c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Simulations\run.jl:130
 [19] run!(sim::Simulation{…}; pickup::Bool)
    @ Oceananigans.Simulations c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Simulations\run.jl:105
 [20] run!(sim::Simulation{…})
    @ Oceananigans.Simulations c:\Users\xinle\Dropbox\MIT\main_oceananigans\Oceananigans.jl\src\Simulations\run.jl:92
 [21] top-level scope
    @ REPL[32]:1
Some type information was truncated. Use `show(err)` to see complete types.

The reason is that GaussianMask is only defined for (x, y, z)

Oceananigans.jl/src/Forcings/relaxation.jl

Lines 149 to 151 in 81dc6c5

	@inline (g::GaussianMask{:x})(x, y, z) = exp(-(x - g.center)^2 / (2 * g.width^2))
	@inline (g::GaussianMask{:y})(x, y, z) = exp(-(y - g.center)^2 / (2 * g.width^2))
	@inline (g::GaussianMask{:z})(x, y, z) = exp(-(z - g.center)^2 / (2 * g.width^2))

But when some topology is Flat one calls

Oceananigans.jl/src/Forcings/relaxation.jl

Lines 90 to 91 in 81dc6c5

	@inline (f::Relaxation)(x₁, x₂, t, field) =
	f.rate * f.mask(x₁, x₂) * (f.target(x₁, x₂, t) - field)

where mask is only expected to take it 2 arguments.

[glwagner]

Nice catch. While x and z seem solvable with new methods, I'm not sure about the y-direction? It seems there is an ambiguity as to which direction is y (eg it could be x, y or y, z). It seems like a bit of work might be required to do this correctly.

[xkykai]

Yeah. I don't have a quick fix right now, but since GaussianMask strictly depends on only one variable we could probably think of ways to use that.

Or maybe pass grid to GaussianMask?

[glwagner]

Or maybe pass grid to GaussianMask?

That would provide a way to fix it. There is also materialize_forcing, which is called inside the model constructor. This gives you the chance to add grid information "under the hood", which relieves the user of needing to pass the grid in. It does make the job slightly more complicated though, since you need to have a constructor that builds a "stub" (eg it has incomplete information), and then a second constructor that adds grid information.