Analysis of simulation results

[parfenyev]

Hello everyone,

Let's say I run a 2D simulation and the output is the velocity fields u and v. At each point in time, I have two arrays, and they characterize the velocity fields at slightly different points in space.

1. Suppose I want to calculate the value of the Reynolds stress component, i.e. uv. Do I need to manually interpolate the velocity fields (in order to multiply them at one point in space) or is there a ready-made solution built into the package?
2. A similar question concerns integration over the entire area of the domain - if I originally had an inhomogeneous grid, is there a ready-made solution for calculating such integrals (e.g., the total kinetic energy)?
3. Finally, quite often it is necessary to carry out analysis in Fourier space - are there any ready-made solutions for calculating, e.g., the energy spectra?

Thank you in advance!

[glwagner]

Do I need to manually interpolate the velocity fields (in order to multiply them at one point in space) or is there a ready-made solution built into the package?

You can use FieldTimeSeries, illustrated here:

Oceananigans.jl/examples/two_dimensional_turbulence.jl

Line 110 in ee0bb44

ω_timeseries = FieldTimeSeries("two_dimensional_turbulence.jld2", "ω")

So this looks like

ut = FieldTimeSeries(filename, "u") # if you called the output "u"
vt = FieldTimeSeries(filename, "v")

n = 1 # time index
u = ut[n]
v = vt[n]

uv = Field(u * v)
compute!(uv)

A similar question concerns integration over the entire area of the domain - if I originally had an inhomogeneous grid, is there a ready-made solution for calculating such integrals (e.g., the total kinetic energy)?

You can use Integral:

U = Field(Integral(u, dims=(1, 2, 3))

Finally, quite often it is necessary to carry out analysis in Fourier space - are there any ready-made solutions for calculating, e.g., the energy spectra?

There aren't tools for spectral analysis in Oceananigans. Perhaps others can chime in on what's out there?

[glwagner]

I also think that this bug should be fixed.

Which bug? They should all be fixed, for sure!

[parfenyev]

Thanks! There is no bug, I overlooked the option boundary_conditions in the FieldTimeSeries. It would probably be good to add an appropriate comment to the boundary conditions section in the docs.

I've got one more question. In my example, I am driving the fluid with an external force:

u_forcing(x, y, z, t) = 0.03*sin(5*y)
v_forcing(x, y, z, t) = -0.03*sin(5*x)

After finishing the calculation, I want to calculate the power that the external force develops. So, I need to calculate the mean value of u_forcing[x,y,t]*u[x,y,t] + v_forcing[x,y,t]*v[x,y,t]. How can this be done effectively? Probavly, the Field should be used, but I don't know how to pass values to it that depend on spatial coordinates.

[glwagner]

It would probably be good to add an appropriate comment to the boundary conditions section in the docs.

I agree, that makes sense to include an example setting boundary conditions for FieldTimeSeries, in addition to Field.

So, I need to calculate the mean value of u_forcing[x,y,t]*u[x,y,t] + v_forcing[x,y,t]*v[x,y,t]. How can this be done effectively?

You can wrap u_forcing and v_forcing in FunctionField (using model.clock to capture the time), and form a new field representing that calculation for output.

Because you're on a staggered grid, you may not want to add those two calculations together immediately. If you're trying to write down a budget, you'll want to first integrate over the domain, and then add the scalars to obtain the total power input. You also need to make sure locations are correct, eg

u_forcing_field = FunctionField(location(u), u_forcing, grid, clock=model.clock)

where u = model.velocities.u.

[parfenyev]

UndefVarError: FunctionField not defined -- Does this option work for version 0.75.3 of Oceananigans.jl?

In my case the force is time-independent, so I end up with:

u_forcing(x, y, z) = 0.03*sin(5*y)
v_forcing(x, y, z) = -0.03*sin(5*x)

u = u_timeseries[5]
v = v_timeseries[5]

u_forcing_field = Oceananigans.Fields.field(location(u), u_forcing, u.grid)
v_forcing_field = Oceananigans.Fields.field(location(v), v_forcing, v.grid)

power = mean(compute!(Field(u*u_forcing_field)))+mean(compute!(Field(v*v_forcing_field)))

[glwagner]

It's not exported at the top level. Try:

using Oceananigans.Fields: FunctionField

But what you've written is just as good (field should build a function behind the scenes).

You only need to pass clock if you want to use the same function for forcing and for diagnostics (ie so you don't have to define it twice).

[loganpknudsen]

Hi, I am trying to implement a time varying Function Field for my diagnostic calculations for the function":

u_pert(x,z,t,p) = p.uₒ*cos(p.fˢ*t)+p.a1*sin(p.fˢ*t)

If I try to write a FunctionField so that it has a time varying component

UPERT = Oceananigans.Fields.FunctionField(location(uw), u_pert, grid, model.clock)

I get the following error in the REPL

ERROR: MethodError: no method matching FunctionField(::Tuple{…}, ::typeof(u_pert), ::RectilinearGrid{…}, ::Clock{…})

Closest candidates are:
  FunctionField(::Tuple, ::Any, ::Any)
   @ Oceananigans ~/.julia/packages/Oceananigans/thcPW/src/Fields/function_field.jl:54

Stacktrace:
 [1] top-level scope
   @ REPL[105]:1
Some type information was truncated. Use `show(err)` to see complete types.

but if I instead remove the time varying component

UPERT = Oceananigans.Fields.FunctionField(location(uw), u_pert, grid)

there is no error until the model starts to run at which point I get the following error:

[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (801.511 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (9.158 seconds).
ERROR: MethodError: no method matching u_pert(::Float64, ::Float64)

Closest candidates are:
  u_pert(::Any, ::Any, ::Any, ::Any)
   @ Main REPL[52]:1

Stacktrace:
  [1] call_func
    @ ~/.julia/packages/Oceananigans/thcPW/src/Fields/function_field.jl:62 [inlined]
  [2] getindex
    @ ~/.julia/packages/Oceananigans/thcPW/src/Fields/function_field.jl:64 [inlined]
  [3] identity3(i::Int64, j::Int64, k::Int64, grid::RectilinearGrid{…}, c::FunctionField{…})
    @ Oceananigans.Operators ~/.julia/packages/Oceananigans/thcPW/src/Operators/interpolation_utils.jl:39
  [4] (::Oceananigans.AbstractOperations.var"#13#14"{…})(γ::Int64)
    @ Oceananigans.AbstractOperations ~/.julia/packages/Oceananigans/thcPW/src/AbstractOperations/multiary_operations.jl:17
  [5] ntuple
    @ ./ntuple.jl:50 [inlined]
  [6] getindex
    @ ~/.julia/packages/Oceananigans/thcPW/src/AbstractOperations/multiary_operations.jl:16 [inlined]
  [7] identity3
    @ ~/.julia/packages/Oceananigans/thcPW/src/Operators/interpolation_utils.jl:39 [inlined]
  [8] +
    @ ~/.julia/packages/Oceananigans/thcPW/src/AbstractOperations/binary_operations.jl:63 [inlined]
  [9] getindex
    @ ~/.julia/packages/Oceananigans/thcPW/src/AbstractOperations/binary_operations.jl:23 [inlined]
 [10] cpu__compute!
    @ ~/.julia/packages/KernelAbstractions/MAxUm/src/macros.jl:287 [inlined]
 [11] cpu__compute!(__ctx__::KernelAbstractions.CompilerMetadata{…}, data::OffsetArrays.OffsetArray{…}, operand::Oceananigans.AbstractOperations.BinaryOperation{…})
    @ Oceananigans.AbstractOperations ./none:0
Some type information was truncated. Use `show(err)` to see complete types.

where it appears it throws an error for not having input for time or parameters. How do I properly include these as inputs for a FunctionField? I can see in the source code there should still be support for it but I am unclear on what I am doing wrong with regards to defining them or calling the version of FunctionField that supports time varying functions with parameters.

[glwagner]

Hmm, I'm not sure where your first signature comes from.

Here's the docstring:

help?> Oceananigans.Fields.FunctionField
  FunctionField{LX, LY, LZ}(func, grid; clock=nothing, parameters=nothing) where {LX, LY, LZ}

  Returns a FunctionField on grid and at location LX, LY, LZ.

  If clock is not specified, then func must be a function with signature func(x, y, z). If clock is specified, func must be a
  function with signature func(x, y, z, t), where t is internally determined from clock.time.

  A FunctionField will return the result of func(x, y, z [, t]) at LX, LY, LZ on grid when indexed at i, j, k.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

  FunctionField{LX, LY, LZ}(func::FunctionField, grid; clock) where {LX, LY, LZ}

  Adds clock to an existing FunctionField and relocates it to (LX, LY, LZ) on grid.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

  FunctionField(L::Tuple, func, grid)

  Returns a stationary FunctionField on grid and at location L = (LX, LY, LZ), where func is callable with signature func(x, y, z).

So you need to use

UPERT = FunctionField{Face, Center, Center}(u_pert, grid, clock=model.clock, parameters=u_pert_parameters)

Is there any way we can help finding the docstring easier? Or interpreting it?

[loganpknudsen]

Ah that makes more sense. I think the main source of confusion was that I though I had to put the location(...) first instead of being able to include it as brackets before the function input. I am going to chop this up to inexperience on my part, but it could be helpful to note that using 'location' calls to a version of the function without time variability and you have to include {LX,LY,LZ} to access time variability(assuming I understand the difference correctly).

[glwagner]

Hm yeah, I don't think we should be using the form that has location as a positional argument. We should always use FunctionField{LX, LY, LZ}.

[loganpknudsen]

I will keep that in mind when writing future code with FunctionField. Thanks again for the help!

[glwagner]

Here's a PR that gets rid of the extra docstrings: #3678

[glwagner]

Also FYI :-D https://www.merriam-webster.com/dictionary/chalk%20up

[loganpknudsen]

Thank you for your help!