Multidimensional Positive Definite Advection Transport Algorithm Research Articles

On numerical broadening of particle-size spectra: a condensational growth study using PyMPDATA 1.0

Abstract. This work discusses the numerical aspects of representing the condensational growth of particles in models of aerosol systems such as atmospheric clouds. It focuses on the Eulerian modelling approach, in which fixed-bin discretisation is used for the probability density function describing the particle-size spectrum. Numerical diffusion is inherent to the employment of the fixed-bin discretisation for solving the arising transport problem (advection equation describing size spectrum evolution). The focus of this work is on a technique for reducing the numerical diffusion in solutions based on the upwind scheme: the multidimensional positive definite advection transport algorithm (MPDATA). Several MPDATA variants are explored including infinite-gauge, non-oscillatory, third-order terms and recursive antidiffusive correction (double-pass donor cell, DPDC) options. Methodologies for handling coordinate transformations associated with both particle-size spectrum coordinate choice and with numerical grid layout choice are expounded. Analysis of the performance of the scheme for different discretisation parameters and different settings of the algorithm is performed using (i) an analytically solvable box-model test case and (ii) the single-column kinematic driver (“KiD”) test case in which the size-spectral advection due to condensation is solved simultaneously with the advection in the vertical spatial coordinate, and in which the supersaturation evolution is coupled with the droplet growth through water mass budget. The box-model problem covers size-spectral dynamics only; no spatial dimension is considered. The single-column test case involves a numerical solution of a two-dimensional advection problem (spectral and spatial dimensions). The discussion presented in the paper covers size-spectral, spatial and temporal convergence as well as computational cost, conservativeness and quantification of the numerical broadening of the particle-size spectrum. The box-model simulations demonstrate that, compared with upwind solutions, even a 10-fold decrease in the spurious numerical spectral broadening can be obtained by an apt choice of the MPDATA variant (maintaining the same spatial and temporal resolution), yet at an increased computational cost. Analyses using the single-column test case reveal that the width of the droplet size spectrum is affected by numerical diffusion pertinent to both spatial and spectral advection. Application of even a single corrective iteration of MPDATA robustly decreases the relative dispersion of the droplet spectrum, roughly by a factor of 2 at the levels of maximal liquid water content. Presented simulations are carried out using PyMPDATA – a new open-source Python implementation of MPDATA based on the Numba just-in-time compilation infrastructure.

A Validated Study of a Modified Shallow Water Model for Strong Cyclonic Motions and Their Structures in a Rotating Tank

A joint theoretical and numerical study was carried out to investigate the fluid dynamical aspect of the motion of a vortex generated in a rotating tank with a sloping bottom. This study aims at understanding the evolution of strong cyclonic motions on a β-plane in the Northern Hemisphere. The strong cyclonic vortices were characterized by four nondimensional parameters which were derived through a scale analysis of the depth variations of fluid. By simplifying the model flow field and the prototype flow field, respectively, through the conservation of potential vorticity, two sets of dynamic similarity conditions are derived. This study proposed a sophisticated modified shallow water model (MSWM) to investigate the flow features of such strong vortices. A detailed numerical calculation adopted by multidimensional positive definite advection transport algorithm (MPDATA) was carried out to validate those effects considered in the MSWM model, including sloping bottom, parabolic free surface deformation, and viscous dissipation. Close agreements were found between the experimental and numerical results, including the streamlines patterns and the vortex trajectory. Comprehensive simulations for strong cyclonic vortices over different sloping bottoms were investigated to understand the impact of planetary β effect on vortex. The results calculated by MSWM demonstrate a variety of flow features of interactions between the primary vortex and induced secondary Rossby wave wakes that were essential and prominent in environmental geophysical flows.

Positive Definite Advection Transport Algorithm for Conservation Law Equations on Nonuniform Irregular Grids

The multidimensional positive definite advection transport algorithm (MPDATA) is an important numerical method for the computation of atmospheric dynamics. MPDATA is second-order accurate, positive definite, conservative, and computationally efficient. However, the method is problematic in which it results in a loss of precision when computing a nonuniform irregular grid. Furthermore, research revealed two reasons for this problem. On the one hand, numerical discretization of boundary derivatives of the finite-volume method is incompatible with nonuniform meshes (or grids); on the other hand, the up-wind scheme of staggered grids is not applicable to the calculation of irregular grids. We overcome these two problems by using the multipoint Taylor expansion method to obtain a boundary derivative numerical approximation scheme that does not depend on the grid structure. Furthermore, combined with the well-balance central-upwind scheme, a positive definite advection scheme for irregular meshes is proposed. Then, the positivity of the new numerical scheme is analyzed. Finally, the result of this study is verified by numerical simulation.

A Multidimensional Positive Definite Remapping Algorithm for Unstructured Meshes

We report on developments of a second-order, conservative, sign-preserving remapping scheme for Arbitrary Lagrangian-Eulerian (ALE) methods operating on unstructured meshes. The remapping uses concepts of the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA).The non-oscillatory infinite gauge option of MPDATA remapping is derived in volume coordinates and is based upon a general and compact edge-based data structure, developed for use within an arbitrary finite volume framework. A Flux Corrected Transport style of limiting ensures monotonicity preservation, while the construction of volume coordinates utilises median dual polygonal finite volume cells.Theoretical developments are supported by numerical testing involving idealised cases with prescribed mesh movement for advection of scalars. The numerical investigations include an asymptotic mesh convergence study and comparisons with both MPDATA and Van Leer MUSCL remapping schemes operating on Cartesian meshes. The results demonstrate that the proposed scheme is suitable for providing accurate ALE remapping for unstructured meshes.

University of Warsaw Lagrangian Cloud Model (UWLCM) 1.0: a modern large-eddy simulation tool for warm cloud modeling with Lagrangian microphysics

Abstract. A new anelastic large-eddy simulation (LES) model with an Eulerian dynamical core and Lagrangian particle-based microphysics is presented. The dynamical core uses the multidimensional positive-definite advection transport algorithm (MPDATA) advection scheme and the generalized conjugate residual pressure solver, whereas the microphysics scheme is based on the super-droplet method. Algorithms for coupling of Lagrangian microphysics with Eulerian dynamics are presented, including spatial and temporal discretizations and a condensation substepping algorithm. The model is free of numerical diffusion in the droplet size spectrum. Activation of droplets is modeled explicitly, making the model less sensitive to local supersaturation maxima than models in which activation is parameterized. Simulations of a drizzling marine stratocumulus give results in agreement with other LES models. It is shown that in the super-droplet method a relatively low number of computational particles is sufficient to obtain correct averaged properties of a cloud, but condensation and collision–coalescence have to be modeled with a time step of the order of 0.1 s. Such short time steps are achieved by substepping, as the model time step is typically around 1 s. Simulations with and without an explicit subgrid-scale turbulence model are compared. Effects of modeling subgrid-scale motion of super-droplets are investigated. The model achieves high computational performance by using graphics processing unit (GPU) accelerators.

Derivative pricing as a transport problem: MPDATA solutions to Black–Scholes-type equations

We discuss applications of the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) to numerically solve partial differential equations arising from stochastic models in quantitative finance. In particular, we develop a framework for solving Black–Scholes-type equations by transforming them into advection–diffusion problems. The equations are then numerically integrated backward in time using an iterative explicit finite-difference approach, in which the Fickian term is represented as an additional advective term. We discuss the correspondence between transport phenomena and financial models, uncovering the possibility of expressing the no-arbitrage principle as a conservation law. We show second-order accuracy in time and space of the embraced numerical scheme. This is done via a convergence analysis comparing MPDATA numerical solutions with classic Black–Scholes analytical formulæ for the valuation of European options. We demonstrate in addition a way of applying MPDATA to solve the free boundary problem (leading to a linear complementarity problem) for the valuation of American options. We finally comment on the potential of MPDATA methods with respect to more complex models typically used in quantitative finance.

Non-oscillatory forward-in-time integrators for viscous incompressible flows past a sphere

A non-oscillatory forward-in-time (NFT) integrator is developed to provide solutions of the Navier–Stokes equations for incompressible flows. Simulations of flows past a sphere are chosen as a benchmark representative of a class of engineering flows past obstacles. The methodology is further extended to moderate Reynolds number, stably stratified flows under gravity, for Froude numbers that typify the characteristic regimes of natural flows past distinct isolated features of topography in weather and climate models. The key elements of the proposed method consist of the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) and a robust non-symmetric Krylov-subspace elliptic solver. The solutions employ a finite volume spatial discretisation on unstructured and hybrid meshes and benefit from a collocated arrangement of all flow variables while being inherently stable. The development includes the implementation of viscous terms with the detached-eddy simulation (DES) approach employed for turbulent flows. Results demonstrate that the proposed methodology enables direct comparisons of the numerical solutions with corresponding laboratory studies of viscous and stratified flows while illustrating accuracy, robustness and flexibility of the NFT schemes. The presented simulations also offer a better insight into stably stratified flows past a sphere.

Performance portable parallel programming of heterogeneous stencils across shared-memory platforms with modern Intel processors

In this work, we take up the challenge of performance portable programming of heterogeneous stencil computations across a wide range of modern shared-memory systems. An important example of such computations is the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA), the second major part of the dynamic core of the EULAG geophysical model. For this aim, we develop a set of parametric optimization techniques and four-step procedure for customization of the MPDATA code. Among these techniques are: islands-of-cores strategy, (3+1)D decomposition, exploiting data parallelism and simultaneous multithreading, data flow synchronization, and vectorization. The proposed adaptation methodology helps us to develop the automatic transformation of the MPDATA code to achieve high sustained scalable performance for all tested ccNUMA platforms with Intel processors of last generations. This means that for a given platform, the sustained performance of the new code is kept at a similar level, independently of the problem size. The highest performance utilization rate of about 41–46% of the theoretical peak, measured for all benchmarks, is provided for any of the two-socket servers based on Skylake-SP (SKL-SP), Broadwell, and Haswell CPU architectures. At the same time, the four-socket server with SKL-SP processors achieves the highest sustained performance of around 1.0–1.1 Tflop/s that corresponds to about 33% of the peak.

Extension of a Non-Oscillatory Forward-in-Time Scheme to Detached Eddy Simulation of Stratified Turbulent Flows With Unstructured Meshes

A detached eddy simulation algorithm is developed based on the non-oscillatory forward-in-time (NFT) multidimensional positive definite advection transport algorithm (MPDATA) scheme for stratified turbulent flows with low or moderate Reynolds numbers on unstructured meshes. NFT MPDATA is extended to viscous turbulent flows for the first time. Viscous diffusion terms and an explicit sub-grid-scale model are included. The working viscosity of the modified Spalart-Allmaras one equation model is transported by an algorithm based on the NFT MPDATA template. The scheme is semi-implicit which solves the linear equation system once per time step. The advantage of the scheme is efficiency and second-order accuracy. Furthermore, an edge-based data structure is adopted such that our scheme is suitable for both structured and unstructured meshes. The linear equation system is solved by an efficient algebraic multigrid (AMG) method. The algorithm is validated using test cases with Reynolds numbers 200 and 5000 over a range of Froude numbers. The flow patterns and dynamic drag coefficients are analysed and compared to numerical and experimental data in literature.

Unleashing the performance of ccNUMA multiprocessor architectures in heterogeneous stencil computations

This paper meets the challenge of harnessing the heterogeneous communication architecture of ccNUMA multiprocessors for heterogeneous stencil computations, an important example of which is the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA). We propose a method for optimization of parallel implementation of heterogeneous stencil computations that is a combination of the islands-of-core strategy and (3{+}1)D decomposition. The method allows a flexible management of the trade-off between computation and communication costs in accordance with features of modern ccNUMA architectures. Its efficiency is demonstrated for the implementation of MPDATA on the SGI UV 2000 and UV 3000 servers, as well as for 2- and 4-socket ccNUMA platforms based on various Intel CPU architectures, including Skylake, Broadwell, and Haswell.

MPDATA: Third-order accuracy for variable flows

This paper extends the multidimensional positive definite advection transport algorithm (MPDATA) to third-order accuracy for temporally and spatially varying flows. This is accomplished by identifying the leading truncation error of the standard second-order MPDATA, performing the Cauchy–Kowalevski procedure to express it in a spatial form and compensating its discrete representation—much in the same way as the standard MPDATA corrects the first-order accurate upwind scheme. The procedure of deriving the spatial form of the truncation error was automated using a computer algebra system. This enables various options in MPDATA to be included straightforwardly in the third-order scheme, thereby minimising the implementation effort in existing code bases. Following the spirit of MPDATA, the error is compensated using the upwind scheme resulting in a sign-preserving algorithm, and the entire scheme can be formulated using only two upwind passes. Established MPDATA enhancements, such as formulation in generalised curvilinear coordinates, the nonoscillatory option or the infinite-gauge variant, carry over to the fully third-order accurate scheme. A manufactured 3D analytic solution is used to verify the theoretical development and its numerical implementation, whereas global tracer-transport benchmarks demonstrate benefits for chemistry-transport models fundamental to air quality monitoring, forecasting and control. A series of explicitly-inviscid implicit large-eddy simulations of a convective boundary layer and explicitly-viscid simulations of a double shear layer illustrate advantages of the fully third-order-accurate MPDATA for fluid dynamics applications.

An unstructured-mesh finite-volume MPDATA for compressible atmospheric dynamics

An advancement of the unstructured-mesh finite-volume MPDATA (Multidimensional Positive Definite Advection Transport Algorithm) is presented that formulates the error-compensative pseudo-velocity of the scheme to rely only on face-normal advective fluxes to the dual cells, in contrast to the full vector employed in previous implementations. This is essentially achieved by expressing the temporal truncation error underlying the pseudo-velocity in a form consistent with the flux-divergence of the governing conservation law. The development is especially important for integrating fluid dynamics equations on non-rectilinear meshes whenever face-normal advective mass fluxes are employed for transport compatible with mass continuity—the latter being essential for flux-form schemes. In particular, the proposed formulation enables large-time-step semi-implicit finite-volume integration of the compressible Euler equations using MPDATA on arbitrary hybrid computational meshes. Furthermore, it facilitates multiple error-compensative iterations of the finite-volume MPDATA and improved overall accuracy. The advancement combines straightforwardly with earlier developments, such as the nonoscillatory option, the infinite-gauge variant, and moving curvilinear meshes. A comprehensive description of the scheme is provided for a hybrid horizontally-unstructured vertically-structured computational mesh for efficient global atmospheric flow modelling. The proposed finite-volume MPDATA is verified using selected 3D global atmospheric benchmark simulations, representative of hydrostatic and non-hydrostatic flow regimes. Besides the added capabilities, the scheme retains fully the efficacy of established finite-volume MPDATA formulations.

Systematic adaptation of stencil‐based 3D MPDATA to GPU architectures

SummaryIn this work, we focus on a systematic adaptation of the stencil‐based multidimensional positive definite advection transport algorithm (MPDATA) to different graphics processing unit (GPU)‐based computing platforms. Another objective of this work is to compare the performance of MPDATA on several platforms, including a multi‐GPU system with two NVIDIA Tesla K80 cards, and single‐card platforms with Tesla K20X, GeForce GTX TITAN, and GeForce GTX 980. The usage of the following optimization methods is proposed to improve the overall performance: (i) reducing the number of operations by the subexpression elimination when implementing 2.5D blocking; (ii) reorganization of boundary conditions for reducing branch instructions; (iii) advanced memory management to increase the coalesced memory access; and (iv) warps rearrangement for optimizing the data access to GPU global memory. The presented methods of the MPDATA adaptation to GPU architectures allow us to efficiently use many graphics processors within a single node by applying peer‐to‐peer data transfers between GPU global memories. We propose an auto‐tuning procedure to compensate architectural differences between the considered platforms. This procedure takes into account algorithm/GPU‐specific parameters. The proposed approach to adaptation of MPDATA to GPU architectures allows us to achieve up to 482.5 Gflop/s for the platform equipped with two NVIDIA K80 GPUs. Copyright © 2016 John Wiley & Sons, Ltd.

Performance modeling of 3D MPDATA simulations on GPU cluster

The goal of this study is to parallelize the multidimensional positive definite advection transport algorithm (MPDATA) across a computational cluster equipped with GPUs. Our approach permits us to provide an extensive overlapping GPU computations and data transfers, both between computational nodes, as well as between the GPU accelerator and CPU host within a node. For this aim, we decompose a computational domain into two unequal parts which correspond to either data dependent or data independent parts. Then, data transfers can be performed simultaneously with computations corresponding to the second part. Our approach allows for achieving 16.372 Tflop/s using 136 GPUs. To estimate the scalability of the proposed approach, a performance model dedicated to MPDATA simulations is developed. We focus on the analysis of computation and communication execution times, as well as the influence of overlapping data transfers and GPU computations, with regard to the number of nodes.

Libmpdata++ 1.0: a library of parallel MPDATA solvers for systems of generalised transport equations

Abstract. This paper accompanies the first release of libmpdata++, a C++ library implementing the multi-dimensional positive-definite advection transport algorithm (MPDATA) on regular structured grid. The library offers basic numerical solvers for systems of generalised transport equations. The solvers are forward-in-time, conservative and non-linearly stable. The libmpdata++ library covers the basic second-order-accurate formulation of MPDATA, its third-order variant, the infinite-gauge option for variable-sign fields and a flux-corrected transport extension to guarantee non-oscillatory solutions. The library is equipped with a non-symmetric variational elliptic solver for implicit evaluation of pressure gradient terms. All solvers offer parallelisation through domain decomposition using shared-memory parallelisation. The paper describes the library programming interface, and serves as a user guide. Supported options are illustrated with benchmarks discussed in the MPDATA literature. Benchmark descriptions include code snippets as well as quantitative representations of simulation results. Examples of applications include homogeneous transport in one, two and three dimensions in Cartesian and spherical domains; a shallow-water system compared with analytical solution (originally derived for a 2-D case); and a buoyant convection problem in an incompressible Boussinesq fluid with interfacial instability. All the examples are implemented out of the library tree. Regardless of the differences in the problem dimensionality, right-hand-side terms, boundary conditions and parallelisation approach, all the examples use the same unmodified library, which is a key goal of libmpdata++ design. The design, based on the principle of separation of concerns, prioritises the user and developer productivity. The libmpdata++ library is implemented in C++, making use of the Blitz++ multi-dimensional array containers, and is released as free/libre and open-source software.

An unstructured-mesh atmospheric model for nonhydrostatic dynamics: Towards optimal mesh resolution

The paper advances the limited-area anelastic model (Smolarkiewicz et al. (2013) [45]) for investigation of nonhydrostatic dynamics in mesoscale atmospheric flows. New developments include the extension to a tetrahedral-based median-dual option for unstructured meshes and a static mesh adaptivity technique using an error indicator based on inherent properties of the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA). The model employs semi-implicit nonoscillatory forward-in-time integrators for soundproof PDEs, built on MPDATA and a robust non-symmetric Krylov-subspace elliptic solver. Finite-volume spatial discretisation adopts an edge-based data structure. Simulations of stratified orographic flows and the associated gravity-wave phenomena in media with uniform and variable dispersive properties verify the advancement and demonstrate the potential of heterogeneous anisotropic discretisation with large variation in spatial resolution for study of complex stratified flows that can be computationally unattainable with regular grids.

Adaptation of MPDATA Heterogeneous Stencil Computation to Intel Xeon Phi Coprocessor

The multidimensional positive definite advection transport algorithm (MPDATA) belongs to the group of nonoscillatory forward-in-time algorithms and performs a sequence of stencil computations. MPDATA is one of the major parts of the dynamic core of the EULAG geophysical model. In this work, we outline an approach to adaptation of the 3D MPDATA algorithm to the Intel MIC architecture. In order to utilize available computing resources, we propose the (3 + 1)D decomposition of MPDATA heterogeneous stencil computations. This approach is based on combination of the loop tiling and fusion techniques. It allows us to ease memory/communication bounds and better exploit the theoretical floating point efficiency of target computing platforms. An important method of improving the efficiency of the (3 + 1)D decomposition is partitioning of available cores/threads into work teams. It permits for reducing inter-cache communication overheads. This method also increases opportunities for the efficient distribution of MPDATA computation onto available resources of the Intel MIC architecture, as well as Intel CPUs. We discuss preliminary performance results obtained on two hybrid platforms, containing two CPUs and Intel Xeon Phi. The top-of-the-line Intel Xeon Phi 7120P gives the best performance results, and executes MPDATA almost 2 times faster than two Intel Xeon E5-2697v2 CPUs.

Adaptation of fluid model EULAG to graphics processing unit architecture

SummaryThe goal of this study is to adapt the multiscale fluid solver EULerian or LAGrangian framewrok (EULAG) to future graphics processing units (GPU) platforms. The EULAG model has the proven record of successful applications, and excellent efficiency and scalability on conventional supercomputer architectures. Currently, the model is being implemented as the new dynamical core of the COSMO weather prediction framework. Within this study, two main modules of EULAG, namely the multidimensional positive definite advection transport algorithm (MPDATA) and the variational generalized conjugate residual, elliptic pressure solver Generalized Conjugate Residual (GCR) are analyzed and optimized. In this paper, a method is proposed, which ensures a comprehensive analysis of the resource consumption including registers, shared, and global memories. This method allows us to identify bottlenecks of the algorithm, including data transfers between host and global memory, global and shared memories, as well as GPU occupancy. We put the emphasis on providing a fixed memory access pattern, padding as well as organizing computation in the MPDATA algorithm. The testing and validation of the new GPU implementation have been carried out based on modeling decaying turbulence of a homogeneous incompressible fluid in a triply‐periodic cube. Simulations performed using the standard version of EULAG and its new GPU implementation give similar solutions. Preliminary results show a promising increase in terms of computational efficiency. Copyright © 2014 John Wiley & Sons, Ltd.

Parallelization of 2D MPDATA EULAG algorithm on hybrid architectures with GPU accelerators

• Parallelization of stencil-based 2D MPDATA algorithm across both GPUs and CPUs. • Minimization of communication between CPU and GPU at the cost of extra computations. • Adaptation of MPDATA to CPUs using space and temporal blocking techniques. • Adaptation of MPDATA to GPUs based on hierarchical decomposition . • Approach to optimization of MPDATA on GPUs using autotuning technique. EULAG (Eulerian/semi-Lagrangian fluid solver) is an established computational model developed for simulating thermo-fluid flows across a wide range of scales and physical scenarios. The dynamic core of EULAG includes the multidimensional positive definite advection transport algorithm (MPDATA) and elliptic solver. In this work we investigate aspects of an optimal parallel version of the 2D MPDATA algorithm on modern hybrid architectures with GPU accelerators, where computations are distributed across both GPU and CPU components. Using the hybrid OpenMP–OpenCL model of parallel programming opens the way to harness the power of CPU–GPU platforms in a portable way. In order to better utilize features of such computing platforms, comprehensive adaptations of MPDATA computations to hybrid architectures are proposed. These adaptations are based on efficient strategies for memory and computing resource management, which allow us to ease memory and communication bounds, and better exploit the theoretical floating point efficiency of CPU–GPU platforms. The main contributions of the paper are: • method for the decomposition of the 2D MPDATA algorithm as a tool to adapt MPDATA computations to hybrid architectures with GPU accelerators by minimizing communication and synchronization between CPU and GPU components at the cost of additional computations; • method for the adaptation of 2D MPDATA computations to multicore CPU platforms, based on space and temporal blocking techniques; • method for the adaptation of the 2D MPDATA algorithm to GPU architectures, based on a hierarchical decomposition strategy across data and computation domains, with support provided by the developed GPU task scheduler allowing for the flexible management of available resources; • approach to the parametric optimization of 2D MPDATA computations on GPUs using the autotuning technique, which allows us to provide a portable implementation methodology across a variety of GPUs. Hybrid platforms tested in this study contain different numbers of CPUs and GPUs – from solutions consisting of a single CPU and a single GPU to the most elaborate configuration containing two CPUs and two GPUs. Processors of different vendors are employed in these systems – both Intel and AMD CPUs, as well as GPUs from NVIDIA and AMD. For all the grid sizes and for all the tested platforms, the hybrid version with computations spread across CPU and GPU components allows us to achieve the highest performance. In particular, for the largest MPDATA grids used in our experiments, the speedups of the hybrid versions over GPU and CPU versions vary from 1.30 to 1.69, and from 1.95 to 2.25, respectively.

On the development and verification of a parametric parallel unstructured-grid finite-volume wind wave model for coupling with ocean circulation models

A parallel, unstructured grid, finite-volume parametric wind-wave model is developed with the intention of coupling to three-dimensional unstructured grid ocean circulation models. The model is derived from a conservation of energy flux formulations. In the model, the shoaling, refraction, and wave dissipation, as well as exchange of current and Stokes drift momentum effects are considered. The unstructured non-overlapping triangle grid cells are adopted in order to fit the curvature of the coastline better. A second-order unwinding scheme in spatial space and a multidimensional positive definite advection transport algorithm (MPDATA) scheme in angle space are used to solve the wave energy equation. The model is parallelized using METIS for domain decomposition with good load balancing across computational nodes. The model has been tested for several idealized situations, and realistic wave simulations for Lake George (near Canberra, Australia) and Gulf of Mexico; the model results agree well with analytical solutions, laboratory observations, and field observations. It is shown that the present model is capable of predicting waves with reasonable accuracy against the observations. The developed model is effective and economical to couple with existing unstructured circulation models.