From my limited knowledge of openfoam and ww3, in my opinion, unless you are very familiar with openfoam, it may be better to use the programming language that you are familiar with, such as C, Fortran, Maltab etc.

Hi Tom,
This sounds like an interesting problem! Ignoring the wave space issue for a moment, I think you need to first cast the problem as a 2D surface using either the finiteArea machinery of 1.6ext, or the thin film model in 2.2.x (see reactingParcelFilmFoam). The wave space issue is a bit more complex and I don't have a good recommendation for an elegant solution. Brute force for a fixed number of wavespace bins would be doable, but ugly. What is the typical wavespace discretization in other codes such as WaveWatch III http://polar.ncep.noaa.gov/waves/wavewatch/ Good luck, Eric 
Hi Tom,
Sounds interesting  I have done waves, but only in 3D. Some clarifications if you don't mind:  you are solving on a 2D surface, right? Is it curved or flat? For the wavenumber space, are we talking 1D, where each wavenumber carries its own Ck, Ctheta pair?  now mane Ck and Ctheta fields will you have? Will it just be N of each, or NxN as you seem to imply? I will call it NxM for the code below... If N, you need PtrList<volScalarField> ks(N); PtrList<volScalarField> thetas(N); You will access them as ks[i] and thetas[i]. If NxM, make a PtrList<PtrList<volScalarField> > bigOne(N); forAll (bigOne, i) { bigOne.set(i, new PtrList<volScalarField>(M)); } and then fill the lot. You will access them as bigOne[i][j]; I have done quite a bit of this sort of thing (OpenFOAM) :) and would be happy to help. I guess you will find my details on the web; if not, please send me a message and we can catch up. Yours, Hrv 
Hello all,
I have been following this thread out of pure curiosity, so allow me to follow up on the equation. The main equation above it indeed in a four dimensional space, where two of the directions are spatial coordinates and the two remaining "directions" is the wave number (or frequency) and the wave propagation direction. Using the above suggested decomposition it is easy to understand the physical system. In each spatial point (x_i, y_i) you have a wave spectrum, where the energy is distributed on a range of frequencies/wave numbers and a range of directions (optimally, the domain will have to be cyclic along the angular axis). The spectrum is represented by the wave action N. This is also why this type of model is called "spectral wave model". Therefore, the solution to the equation is not a deterministic description of the behaviour of the free surface as e.g. obtained with interFoam/waveFoam, but rather a statistical representation in a given point. E.g. in a simple 1D spatial domain, you would have the JONSWAP spectrum (http://www.google.nl/url?sa=i&rct=j&...78322546123846), since the wave can only propagate in one direction (or rather positive and negative directions). The propagation speeds (C) for x and y are simple propagation speeds typically computed from the linear dispersion relation including the effect of wavecurrent interaction. The C_theta is a propagation speed related to refraction, i.e. advection of wave action from one propagation direction to another due to changes in the bathymetry. I am not quite sure with respect to the last velocity. The source terms on the right hand side include effects from bottom friction, wind drag, dissipation due to wave breaking and nonlinear transfer functions from one frequency to another. The latter involves either three or four frequencies/wave numbers (triad and quad interactions of the wave components). Finally, to answer the question of the layout of the computational domain. The well know commercial and opensource tools available (Wave Watch III, SWAN and MIKE SW) all solve these equations on flat meshes, i.e. local conditions around a harbour, or on spherical meshes for the global prediction of wave conditions. This type of models are then typically coupled to large scale hydrodynamic models/oceanographics models, such that predictions of storm surge, wave induced circulations, etc can be simulated. All of these parameters have importance for a wide variety of topics: swimmers safely, coastal protection, disaster management (should we or should we not evacuate this area), navigation, etc. Kind regards Niels 
Niels, nice summary!
Quote:

Hi all, In addition to the summary of Niels, I would like to add the following: Presently I’m working on solver that solves the above mentioned equation, neglecting the variation in frequency space from the above shown wave action equation (so that results in the dimensions time, x, y and direction), similar to Holthuijsen et al (1989). The implementation of that wave action equation was possible using readily available components from OpenFOAM and assigning the directional space to the zaxis. Results are unfortunately not yet ready for publication, that will take some more time, but it seems to be working. The next step is to extend the above mentioned solver with the frequency space. For this, a “trick” is needed to allow for the extra dimension within OpenFOAM. One way is, as above mentioned by Mr. Jasak, to create additional fields for the different frequencies to be considered (or even extra fields for each frequency and direction component if you use the FaMesh approach from Mr. Patterson). This would subsequently also require a connectivety table that can be used for solving energy transport from one frequency to another (either by means of the propagation or by means of source terms). There are also other “work arounds” to obtain the same result: either by splitting up the action equation (as stated in the first post), using multiple meshes, or using a VectorN approach. I will be working on the next step (including the frequency space) in a few months time, using approaches from the spectral wave models WaveWatch and SWAN. In this next step I will use specific schemes for discretisation of the wave action equation that we are presently investigating as part of our long term research and development of SWAN. I’ll try and keep you posted on the progress/results. Regards, Alfons Smale Holthuijsen, L.H., Booij, N., Herbers, T.H.C., 1989. A prediction model for stationary shortcrested waves in shallow water with ambient currents. Coast. Eng. 13, 23–54. 
Thank you all for your interest in this problem and the suggestions.
Eric, Thanks for the suggestions on finiteArea and the thin film model. I'll look into them. As far as wavespace discretization schemes, here's a summary of what I've seen: WaveWatch III has both a simple first order upwind scheme as well as a third order / limited scheme available (which they call ULTIMATE QUICKEST). Explicit methods are used to solve the equation. SWAN uses a hybrid central difference  upwind scheme. In this case the equation is solved with an implicit method. (In both WaveWatch III and SWAN the same methods are applied to both wavenumber magnitude and direction space.) There are also a couple of unstructured implementations of SWAN; one using a finite element method (Hsu, Ou and Liau, 'Hindcasting nearshore wind waves using a FEM code for SWAN', Coastal Engineering 52 (2005) 177195) and another using a finite volume method (Qi, Chen, Beardsley, Perrie, Cowles and Lai, 'An unstructuredgridfinitevolume surface wave model (FVCOMSWAVE): Implementation, validations and applications', Ocean Modelling 28 (2009) 153166). Both use the Flux Corrected Transport scheme in the wavenumber magnitude space and a central differencing / Crank Nicholson scheme in wavenumber direction space. In all cases, a logarithmic distribution discretization is used for the wavenumber magnitude and a uniform distribution in wavenumber direction. My initial impression is in line with yours; a brute force method should be doable. My plan is to start with that and work towards something more elegant. Hrvoje, yes, it is a 2D surface. Both WaveWatch III and SWAN solve a spherical coordinate form of the equations allowing for curved surfaces. My plan is to solve the Cartesian form and just consider flat surfaces. "For the wavenumber space, are we talking 1D, where each wavenumber carries its own Ck, Ctheta pair?" Yes, each point in wavespace (i.e., each wavenumber magnitude and direction pair) has a distinct Ck, Ctheta pair. However, Ck and Ctheta are also functions of your position in geographical space since they are a function of the surface current speeds and ocean depth. This is also true for C_x and C_y. You can store k and theta with 1D arrays since the discretization does not vary with geographical space (as you suggest in the code contained in your response). Also note that all of the C's are only functions of the surface currents, ocean depth and wavespace position. So they need only be calculated up front and stored for use through the solution procedure. I believe I need to store these things with 3D arrays; one for geographical space and two for wavespace. For ks and thetas your suggestion works: PtrList<volScalarField> ks(N); PtrList<volScalarField> thetas(M); But, for the C's and the wave action I need C(I,N,M) where I identifies the position in geographical space. Is there a way to generate an array that can be accessed bigOne[i][n][m], where i is the mesh index, n is the wavenumber magnitude index and m is the wavenumber direction index? If so, for each position in wavespace [n][m], I can copy bigone[i][n][m] into smallOne[i] to solve the second step equations (as outlined in my first post). Your offer to help is greatly appreciated. I will in all likelihood take you up on it as I move forward. Niels and Alfons, thanks for your posts. I am interested in following your progress Alfons. I didn't know there was a version of the equations that included just the wave direction and not the frequency. 
All times are GMT 4. The time now is 21:15. 