Edit on GitHub

How to...

WaterWaves1D.jl is meant to be versatile, and integrating new blocks to the package is easy. If you ever do so, please do not hesitate to contact the developers to either get help, or report on your advances.

build your model

Let us add the linear (Airy) water waves model whose equations are (using the notations introduced here)

\[ \left\{\begin{array}{l} ∂_tη+\frac{\tanh(\sqrt{μ} D)}{ν\sqrt{μ} D}∂_xv=0,\\[1ex] ∂_tv+∂_xη=0, \end{array}\right.\]

where we use the notation $F(D)$ for the action of pointwise multiplying by the function $F$ in the Fourier space.

In a dedicated file we write

export Airy
mutable struct Airy <: AbstractModel
  label   :: String
  f!      :: Function
  mapto   :: Function
  mapfro  :: Function

  # We build our model here.

end

Our purpose is to provide

  • label, a string used in subsequent informational messages, plots, etc.
  • f!, a function to be called in explicit time-integration solvers such as Euler or RK4 (one may provide other functions to be used with other solvers such as EulerSymp)
  • mapto, a function which from a couple (η,v) of type InitialData provides the raw data matrix on which computations are to be executed
  • mapfro, a function which from raw data returns (η,v,x) where
    • η is the values of surface deformation at collocation points x;
    • v is the derivative of the trace of the velocity potential at x.

To this aim, in place of the commented line, we write

function Airy(param::NamedTuple; # param is a NamedTuple containing all necessary parameters
  label = "linear (Airy)"  # using a keyword argument allows the user to supersede the default label.
   )

  # Set up
  μ = param.μ
  if !in(:ν,keys(param)) # set default ν if it is not provided
    ν = min(1,1/√μ)
  else
    ν = param.ν
  end
  # Collocation points and Fourier modes
  m = Mesh(param)
  x, k = m.x, m.k
  # Fourier multipliers
  ∂ₓ	 = 1im * k            # Differentiation
  ∂ₓF₁ = 1im * tanh.(√μ*k)/(√μ*ν)
  # Pre-allocation
  fftη = zeros(Complex{Float64}, m.N)
  fftv = zeros(Complex{Float64}, m.N)

  # Evolution equations are ∂t U = f(U)
  function f!(U)
    fftη .= U[:,1]
    fftv .= U[:,2]

    U[:,1] .= -∂ₓF₁.*fftv
    U[:,2] .= -∂ₓ.*fftη
  end

  # Build raw data from physical data (discrete Fourier transform)
  function mapto(data::InitialData)
    U = [fft(data.η(x)) fft(data.v(x))]
  end

  # Return physical data `(η,v,x)` from raw data
  function mapfro(U)
    real(ifft(U[:,1])),real(ifft(U[:,2])),x
  end

  new( label, f!, mapto, mapfro )
end

A useful tool used in the above was the function Mesh. It takes as argument a NamedTuple containing typically a number of points/modes and the (half-)size of the domain and provides the vector of collocations points and Fourier wavenumbers (see this discussion).

The Airy model can now be built as follows

using WaterWaves1D
# include your file
model = Airy((μ=1,L=2π,N=2^8))

build your initial data

The simplest way to build an initial data is to use the function Init, which takes as argument either

  • a function η and a function v (in this order);
  • an array of collocation points and two vectors representing η(x) and v(x) (in this order);
  • a mesh (generated with Mesh) and two vectors representing η(mesh.x) and v(mesh.x) (in this order).

Some relevant initial data (e.g. travelling waves) are built-in; see the library section. You can build your own in the following lines.

struct Heap <: InitialData
	η
	v
	label :: String
	info  :: String

	function Heap(L)
		η=x->exp.(-(L*x).^2)
    v=x->zero(x)
		init = Init(η,v)
		label = "Heap"
		info = "Heap of water, with length L=$L."

		new( init.η,init.v,label,info  )
	end
end

The corresponding initial data can now be built as follows

using WaterWaves1D
# include your file
init = Heap(1)

build your time solver

As an example, let us review how the explicit Euler solver, Euler, is built.

In a dedicated file we write

struct Euler <: TimeSolver
    U1 :: Array
    label :: String

    function Euler( U :: Array; realdata=nothing )
        U1 = copy(U)
        if realdata==true
            U1 = real.(U1);
        end
        if realdata==false
            U1 = complex.(U1);
        end
        new( U1, "Euler" )
    end
end

Here, Euler.U1 is a pre-allocated vector which can be used to speed-up calculations, and Euler-label is the string "Euler", used for future references. The optional keyword argument realdata allows to specify the type of data which the solver will take as arguments: either complex or real vectors.

We shall now add one method to the function step!, performing the explicit Euler step: that is replacing a vector U with U+dt*f(U) where f is provided by the model at stake.

export step!
function step!(solver :: Euler,
                model :: AbstractModel,
                U  ,
                dt )

    solver.U1 .= U        # allocate U to U1
    model.f!( solver.U1 ) # model.f!(U) replaces its argument U with f(U)
    U .+= dt .* solver.U1 # update U
end

access to and manage your data

Once an initial-value problem problem has been solved (i.e. numerically integrated), the raw data is stored in problem.data.U, which is an array whose elements correspond (in chronological order) to values at the different computed times, problem.times.ts. 

param  = ( c = 1.1, ϵ = 1, μ = 1, N  = 2^8, L = 16, T = 1, dt = 0.1)
init = SolitaryWhitham(param)
model = Airy(param)
problem = Problem( model, init, param )
solve!(problem; verbose=false)

problem.data.U
11-element Vector{Matrix{ComplexF64}}:
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -7.002608952054088 + 2.5479481718517273e-9im -6.762869498490983 + 2.460715962160613e-9im; … ; 5.494478771973449 + 3.998408458723641e-9im 5.275526526412351 + 3.839074433673005e-9im; -7.002608952054088 - 2.5479481290623095e-9im -6.762869498490982 - 2.4607158783645465e-9im]
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -7.001276215137731 + 0.13109975548514471im -6.7615823887965 + 0.1374871846894821im; … ; 5.490447810332899 + 0.19709016866706217im 5.2716561965685065 + 0.2157149090370213im; -7.001276215137731 - 0.13109975548514466im -6.761582388796499 - 0.137487184689482im]
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -6.99727851178253 + 0.2621496065989485im -6.757721549743839 + 0.27492203378319585im; … ; 5.4783608404769994 + 0.3938911474404455im 5.260050886428859 + 0.43111330068426723im; -6.99727851178253 - 0.2621496065989485im -6.757721549743838 - 0.2749220337831958im]
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -6.990617363674671 + 0.3930996730616719im -6.751288450923144 + 0.4122522365279363im; … ; 5.458235597339906 + 0.5901141787536514im 5.24072762419929 + 0.6458791296583019im; -6.990617363674671 - 0.3930996730616719im -6.751288450923143 - 0.41225223652793624im]
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -6.9812953063142 + 0.5239001100276662im -6.742285541029758 + 0.5494255195423787im; … ; 5.430101610229897 + 0.7854713490508455im 5.213714762459864 + 0.8596972749813389im; -6.9812953063142 - 0.5239001100276662im -6.742285541029757 - 0.5494255195423786im]
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -6.969315888049911 + 0.6545011296062687im -6.730716246932153 + 0.6863896691751392im; … ; 5.394000159501703 + 0.979676015234227im 5.179051936563736 + 1.0722540061901034im; -6.969315888049911 - 0.6545011296062687im -6.730716246932152 - 0.6863896691751392im]
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -6.954683668728705 + 0.7848530198130697im -6.716584972367529 + 0.8230925513793843im; … ; 5.349984215986777 + 1.1724432252488752im 5.136790006481111 + 1.2832374436653486im; -6.954683668728705 - 0.7848530198130697im -6.716584972367528 - 0.8230925513793843im]
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -6.937404217959926 + 0.914906163492287im -6.699897096265571 + 0.9594821315571391im; … ; 5.298118363270418 + 1.3634901361863676im 5.08699098217357 + 1.4923380162456727im; -6.937404217959926 - 0.914906163492287im -6.69989709626557 - 0.9594821315571391im]
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -6.9174841129953455 + 1.0446110572030451im -6.680658970701013 + 1.0955064943657415im; … ; 5.238478702929748 + 1.552536429293696im 5.029727932608275 + 1.6992489154541872im; -6.9174841129953455 - 1.0446110572030451im -6.680658970701012 - 1.0955064943657415im]
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -6.8949309362256 + 1.1739183300623695im -6.658877918475787 + 1.2311138634789036im; … ; 5.171152742871617 + 1.7393047212785473im 4.96508487854555 + 1.9036665456715618im; -6.8949309362256 - 1.1739183300623695im -6.6588779184757865 - 1.2311138634789036im]
 [7.666279752905108 + 0.0im 7.418980406037189 + 0.0im; -6.869753272294022 + 1.3027787625377236im -6.634562230331669 + 1.3662526212948567im; … ; 5.096239268934254 + 1.9235209713074486im 4.893156669257143 + 2.1052909695949173im; -6.869753272294022 - 1.3027787625377236im -6.634562230331668 - 1.3662526212948567im]

Physical data (say at final computed time) can be reconstructed using the model as follows: 

η,v,x = problem.model.mapfro(last(problem.data.U))
([1.3452999706925484e-5, 1.2648564822109742e-5, 1.1943588073548561e-5, 1.1332528170621137e-5, 1.081058232824117e-5, 1.0373648571218586e-5, 1.0018293550059787e-5, 9.741725618517272e-6, 9.541772960947048e-6, 9.4168666037453e-6  …  2.882170424055872e-5, 2.6515672495059706e-5, 2.4418100230207823e-5, 2.2512500854929285e-5, 2.0783896173305616e-5, 1.9218698806965084e-5, 1.7804605523853018e-5, 1.653050064017239e-5, 1.5386368739392364e-5, 1.436321602723764e-5], [1.3452888219009362e-5, 1.2797168146588822e-5, 1.2242070555210094e-5, 1.1783230469991635e-5, 1.1417039629257886e-5, 1.114061811510747e-5, 1.0951791690416712e-5, 1.0849074692023675e-5, 1.0831658321108772e-5, 1.0899404251002442e-5  …  2.713429311568087e-5, 2.503262306819895e-5, 2.3127772713633487e-5, 2.1404766900439895e-5, 1.9850059768646705e-5, 1.8451428328175937e-5, 1.7197876434071196e-5, 1.6079548373168495e-5, 1.508765143203028e-5, 1.4214386786576966e-5], [-16.0, -15.875, -15.75, -15.625, -15.5, -15.375, -15.25, -15.125, -15.0, -14.875  …  14.75, 14.875, 15.0, 15.125, 15.25, 15.375, 15.5, 15.625, 15.75, 15.875])

where η and v are respectively the values of the surface deformation and of the derivative of the trace of the velocity potential at the surface, at collocation points x.

This procedure is carried out by the function solution, which allows in addition to perform some interpolations (making use of the otherwise helpful function interpolate).

η,v,x = solution(problem)
([1.3452999706925484e-5, 1.2648564822109742e-5, 1.1943588073548561e-5, 1.1332528170621137e-5, 1.081058232824117e-5, 1.0373648571218586e-5, 1.0018293550059787e-5, 9.741725618517272e-6, 9.541772960947048e-6, 9.4168666037453e-6  …  2.882170424055872e-5, 2.6515672495059706e-5, 2.4418100230207823e-5, 2.2512500854929285e-5, 2.0783896173305616e-5, 1.9218698806965084e-5, 1.7804605523853018e-5, 1.653050064017239e-5, 1.5386368739392364e-5, 1.436321602723764e-5], [1.3452888219009362e-5, 1.2797168146588822e-5, 1.2242070555210094e-5, 1.1783230469991635e-5, 1.1417039629257886e-5, 1.114061811510747e-5, 1.0951791690416712e-5, 1.0849074692023675e-5, 1.0831658321108772e-5, 1.0899404251002442e-5  …  2.713429311568087e-5, 2.503262306819895e-5, 2.3127772713633487e-5, 2.1404766900439895e-5, 1.9850059768646705e-5, 1.8451428328175937e-5, 1.7197876434071196e-5, 1.6079548373168495e-5, 1.508765143203028e-5, 1.4214386786576966e-5], [-16.0, -15.875, -15.75, -15.625, -15.5, -15.375, -15.25, -15.125, -15.0, -14.875  …  14.75, 14.875, 15.0, 15.125, 15.25, 15.375, 15.5, 15.625, 15.75, 15.875], 1.0)

Using (η,v,x) one can compute other quantities such as the mass, momentum, energy; for instance for the purpose of testing how well these quantities are numerically perserved (when the quantities are first integrals of the considered model). Built-in functions mass, momentum, energy (and massdiff, momentumdiff, energydiff) compute such quantities (or their variation).

plot your data

Once (η,v,x) is obtained as above, producing plots is as simple as

using Plots
plot(x, [ η v ], label = ["η" "v"])

One can also plot the amplitude of discrete Fourier coefficients (in semi-log scale) as follows:

using FFTW
k=fftshift(Mesh(x).k);fftη=fftshift(fft(η));fftv=fftshift(fft(v));
indices = (fftη .!=0) .& (fftv .!=0 )
plot(k[indices], [ abs.(fftη)[indices] abs.(fftv)[indices] ], yscale=:log10)

Built-in plot recipes provide convenient ways of producing such plots, and animations.

plot(problem;var=[:surface,:velocity,:fourier])

save and load your data

Thanks to the package HDF5.jl it is possible to save (raw) data to a local file, and then load them for future analyses.

Built-in functions ease the process. Here is a typical example.

Build and solve a problem.

param  = ( c = 1.1, ϵ = 1, μ = 1, N  = 2^8, L = 16, T = 1, dt = 0.1)
init = SolitaryWhitham(param)
model = Airy(param)
problem = Problem( model, init, param )
solve!(problem; verbose=false)
nothing
[ Info: Converged : relative error 3.517149510517202e-15 in 4 steps

Save the initial data.

x = Mesh(param).x
dump("file_name", x, init)

Save the time-integrated (raw) data.

dump("file_name", problem) # or dump("file_name", problem.data)

Reconstruct the problem, without solving it.

loaded_init = load_init("file_name") # load initial data
new_problem = Problem( model, loaded_init, param ) # re-build problem
load_data!("file_name", new_problem) # incorporate time-integrated data

problem.data == new_problem.data
true
Note

Initial data are functions. The dump function saves values at some collocation points. Hence the loaded initial data typically differs from the original one by machine epsilon rounding errors.

Note

It is not possible to save models, solvers, or problems. Hence the user needs to store separately the parameters and information required to build them.