Accelerating Large‐Eddy Simulations of Clouds With Tensor Processing Units - J Adv Model Earth Syst - 2023 - Chammas - Accelerating Large‐Eddy Simulations of Clouds With Tensor Processing Units.pdf

Introduction

Scientific progress accelerates when it is possible to cycle rapidly through the knowledge discovery loop: design

and conduct experiments, learn from the experiments, and design and conduct new experiments to test and

refine models and hypotheses with the information obtained from them (National Academies of Sciences, Engi

neering, and Medicine,

2022

). In the computational sciences, experiments are conducted numerically, and the

ability to cycle through the knowledge discovery loop has advanced hand-in-hand with the evolution of computer

hardware. The atmospheric sciences represent a prime example of advances in computer hardware enabling

and accelerating scientific progress. The first experiments with two-dimensional atmosphere models (Charney

et al.,

1950

) and, soon thereafter, with quasigeostrophic two-layer models (Phillips,

1954

1956

) only allowed

simulations that were slower than or comparable with the real-time evolution of the atmosphere. The first experi

ments using general circulation models similarly pushed the envelope of what was computationally feasible at the

time (Manabe et al.,

1965

; Smagorinsky,

1963

; Smagorinsky et al.,

1965

). Once such simulations of atmospheric

flows, albeit at coarse resolution, became routine and rapidly executable, systematic exploration and experimen

tation followed, enabling rapid progress in our understanding of the atmosphere's general circulation, from its

dependence on planetary characteristics such as planetary radius and rotation rate (Williams,

1988a

1988b

), over

the nature of atmospheric turbulence (Held,

1999

; Held & Larichev,

1996

; Rhines,

1975

1979

; Schneider,

2006

;

Abstract

Clouds, especially low clouds, are crucial for regulating Earth's energy balance and mediating

the response of the climate system to changes in greenhouse gas concentrations. Despite their importance

for climate, they remain relatively poorly understood and are inaccurately represented in climate models. A

principal reason is that the high computational expense of simulating them with large-eddy simulations (LES)

has inhibited broad and systematic numerical experimentation and the generation of large data sets for training

parametrization schemes for climate models. Here we demonstrate LES of low clouds on tensor processing

units (TPUs), application-specific integrated circuits that were originally developed for machine learning

applications. We show that TPUs in conjunction with tailored software implementations can be used to simulate

computationally challenging stratocumulus clouds in conditions observed during the Dynamics and Chemistry

of Marine Stratocumulus (DYCOMS) field study. The TPU-based LES code successfully reproduces clouds

during DYCOMS and opens up the large computational resources available on TPUs to cloud simulations.

The code enables unprecedented weak and strong scaling of LES, making it possible, for example, to simulate

stratocumulus with 10× speedup over real-time evolution in domains with a 34.7 km × 53.8 km horizontal

cross section. The results open up new avenues for computational experiments and for substantially enlarging

the

sample of LES available to train parameterizations of low clouds.

Plain Language Summary

The study of clouds has been impeded by, among other factors,

limitations in our ability to simulate them rapidly and on sufficiently large domains. In particular,

computational limitations in simulating low clouds are among the reasons for the difficulties of representing

them accurately in climate models; this is one of the dominant uncertainties in climate predictions. This

paper demonstrates how the large computing power available on tensor processing units (TPUs) (integrated

circuits originally designed for machine learning applications) can be harnessed for simulating low clouds.

We demonstrate the largest simulations of low clouds to date, with hundreds of billions of variables, and we

document their fidelity to aircraft observations. The results open up the large computational resources available

on TPUs, hitherto primarily used for machine learning, to the study of clouds in the climate system.

CHAMMAS ET AL.

in Modeling Earth Systems published

by Wiley Periodicals LLC on behalf of

American Geophysical Union.

This is an open access article under

the terms of the

Creative Commons

Attribution

License, which permits use,

distribution and reproduction in any

medium, provided the original work is

properly cited.

Accelerating Large-Eddy Simulations of Clouds With Tensor

Processing Units

Sheide Chammas

, Qing Wang

, Tapio Schneider

1,2

, Matthias Ihme

1,3

, Yi-fan Chen

, and

John Anderson

Google LLC, Mountain View, CA, USA,

California Institute of Technology, Pasadena, CA, USA,

Stanford University,

Stanford, CA, USA

Key Points:

•

We introduce a large-eddy simulation

(LES) framework that runs on tensor

processing units (TPUs, accelerators

designed for machine learning)

•

The fidelity of the LES is established

by reproducing aircraft observations

of nocturnal stratocumulus clouds

over the Pacific

•

The LES exhibit unprecedented

scalability on TPUs, enabling the

large-scale generation of training data

for cloud parameterizations

Correspondence to:

S. Chammas,

sheide@google.com

Citation:

Chammas, S., Wang, Q., Schneider,

T., Ihme, M., Chen, Y.-f., &

Anderson, J. (2023). Accelerating

large-eddy simulations of clouds with

tensor processing units.

Journal of

Advances in Modeling Earth Systems

, e2023MS003619.

https://doi.

org/10.1029/2023MS003619

Received 23 JAN 2023

Accepted 13 SEP 2023

Author Contributions:

Conceptualization:

Tapio Schneider,

John Anderson

Formal analysis:

Sheide Chammas,

Tapio Schneider, Yi-fan Chen

Investigation:

Sheide Chammas, Tapio

Schneider, Matthias Ihme

Methodology:

Matthias Ihme

Project Administration:

Yi-fan Chen,

John Anderson

Software:

Sheide Chammas, Qing Wang,

Yi-fan Chen

Supervision:

Tapio Schneider, John

Anderson

Validation:

Sheide Chammas, Qing

Wang, Tapio Schneider, Matthias Ihme,

John Anderson

10.1029/2023MS003619

RESEARCH ARTICLE

1 of 18

Journal of Advances in Modeling Earth Systems

CHAMMAS ET AL.

10.1029/2023MS003619

2 of 18

Schneider & Walker,

2006

), to elucidating the hydrologic cycle (Allen & Ingram,

2002

; Chou & Neelin,

2004

;

Held & Soden,

2006

; Manabe & Wetherald,

1975

; O’Gorman & Schneider,

2008

; Rind et al.,

1992

; Schneider

et al.,

2010

). Similarly, our understanding of deep convective clouds advanced substantially once deep-convection

resolving simulations in limited areas became routinely feasible (T. Cronin,

2014

; T. W. Cronin et al.,

2015

;

Held et al.,

1993

; Tompkins & Craig,

1998

; Wing et al.,

2018

). In contrast, our understanding of the dynam

ics of low clouds is in its infancy. We do not have quantitative theories of their response to climate change

(Bretherton,

2015

), and shortcomings in their representation in climate models have long dominated uncer

tainties in climate projections (Bony & Dufresne,

2005

; Brient & Schneider,

2016

; Brient et al.,

2016

; Cess

et al.,

1990

1996

; Dufresne & Bony,

2008

; Vial et al.,

2013

; Webb et al.,

2006

2013

; Zelinka et al.,

2017

Numerical experiments have been limited to studies that have explored a few dozen canonical situations,

mostly in the tropics (Blossey et al.,

2013

2016

; Caldwell & Bretherton,

2009

; Rauber et al.,

2007

; Sandu &

Stevens,

2011

; Schalkwijk et al.,

2015

; Siebesma et al.,

2003

; B. Stevens et al.,

2005

; Tan et al.,

2016

2017

;

Zhang et al.,

2012

2013

). Broader exploration has been limited by the computational expense necessary to

resolve the meter-scale dynamics of low clouds in large-eddy simulations (LES).

Here we take the next step in the co-evolution of science and computing hardware by demonstrating LES of low

clouds on tensor processing units (TPUs). TPUs are application-specific integrated circuits, originally developed for

machine learning applications, which are dominated by dense vector and matrix computations (Jouppi et al.,

2017

The current TPU architecture integrates 4,096 chips into a so-called TPU Pod, which achieves 1.1 exaflops in aggre

gate at half precision. TPUs are publicly available for cloud computing and can be leveraged for fluids simulations

(Wang et al.,

2022

) and other scientific computing tasks (Belletti et al.,

2019

; Lu et al.,

2020

; Pederson et al.,

2022

with remarkable computational throughput and scalability. Large, high-bandwidth memory and fast chip-to-chip

interconnects (currently 1.1 PB/s) contribute to the performance of TPUs and alleviate bottlenecks that computa

tional fluid dynamics (CFD) applications typically face on accelerator platforms (Balaji,

2021

). However, the native

half- or single-precision arithmetic of TPUs can also create challenges in CFD applications (Wang et al.,

2022

The objective of this study is to evaluate the throughput and scalability achievable on TPUs in simulations of

subtropical stratocumulus clouds under conditions encountered during the Dynamics and Chemistry of Marine Stra

tocumulus (DYCOMS) field study (B. Stevens et al.,

2005

). Stratocumulus clouds are a particularly good testbed

for low-cloud simulations for two reasons: First, they are the most frequent cloud type on Earth, covering about 20%

of tropical oceans, with an outsize impact on Earth's energy balance (Wood,

2012

). Reductions or increases in the

area they cover by a mere 4% can have an impact on Earth's surface temperature comparable to doubling or halving

atmospheric carbon dioxide concentrations (Randall et al.,

1984

). Second, they are notoriously difficult to simulate,

even in LES, because key processes responsible for their maintenance, such as turbulent entrainment of air across the

often sharp temperature inversions at their tops, occur on scales of meters (Mellado,

2016

). The resulting numerical

challenges lead to large differences among various LES codes owing to differences in the numerical discretizations

(B. Stevens et al.,

2005

; Pressel et al.,

2017

). For example, weighted essentially non-oscillatory (WENO) advection

schemes at resolutions of O(10 m) lead to more faithful simulations—relative to field measurements—than centered

difference advection schemes at resolutions of O(1 m) (Schneider et al.,

2019

, their Supporting Figure 3). These two

reasons make progress in simulating subtropical stratocumulus both important and challenging.

This paper is structured as follows. Section

describes the governing equations, numerical methods, and

TPU-specific implementation decisions in our LES code. Section

presents a dry buoyant bubble and a density

current as validation examples of the code. Section

presents the DYCOMS simulations, including comparisons

with field data and a scaling analysis of the simulations. Section

summarizes the conclusions and new opportu

nities afforded by this TPU-enabled cloud-simulation capability.

Model Formulation, Numerics, and TPU Implementation

2.1.

Governing Equations

Our LES simulates the anelastic equations for moist air, understood to be an ideal admixture of dry air, water vapor,

and any condensed water that is suspended in and moves with the air. Precipitating condensate (e.g., rain and snow)

is not considered part of the working fluid, and the suspended constituents of the moist air are taken to be in local

thermodynamic equilibrium. By Gibbs' phase rule, then, a complete thermodynamic description of this system

with two components (dry air and water) and three phases (water vapor, liquid water, ice) requires specification of

two

thermodynamic variables, in addition to the density

and pressure

of the moist air. We choose the total water

Visualization:

Sheide Chammas, Qing

Wang, Yi-fan Chen

Writing – original draft:

Sheide

Chammas, Qing Wang, Tapio Schneider,

Yi-fan Chen

Writing – review & editing:

Sheide

Chammas, Qing Wang, Tapio Schneider,

Matthias Ihme, Yi-fan Chen

19422466, 2023, 10, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023MS003619, Wiley Online Library on [17/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Advances in Modeling Earth Systems

CHAMMAS ET AL.

10.1029/2023MS003619

3 of 18

specific humidity

(total mass of water per unit mass of moist air) and liquid-

ice potential temperature

(Tripoli & Cotton,

1981

). This choice of thermody

namic variables is advantageous because both the total specific humidity

and

(approximately) the liquid-ice potential temperature

are materially conserved

even in the presence of reversible phase transitions of water. The temperature

and specific humidities

and

of cloud liquid and ice can be computed from

the other thermodynamic variables.

The anelastic approximation eliminates physically insignificant acoustic

waves by linearizing the density

(

) =

(

) +

′(

) and pres

sure

(

) =

(

) +

′(

) around a dry reference state with density

(

) and hydrostatic pressure

(

), which depend only on altitude

. Here,

reference state variables are indicated by a subscript 0, and perturbation vari

ables by primes. Perturbation variables are retained only where they affect

accelerations. The reference density and pressure depend only on the vertical

coordinate

and are in hydrostatic balance,

휕휕휕휕

′

푧푧

)

휕휕푧푧

=−

휌휌

(

푧푧

)

푔푔푔

(1)

For energetic consistency, the reference state needs to be adiabatic, that is,

the reference potential temperature

needs to be constant (Bannon,

1996

;

Pauluis,

2008

). Therefore,

푇푇

′

휃휃

(

푝푝

)

푅푅

푑푑

∕

푐푐

푝푝푑푑

(2)

푝푝

(

1−

푔푔푔푔

푐푐

푝푝푝푝

휃휃

)

푐푐

푝푝푝푝

∕

푅푅

푝푝

(3)

휌휌

′

푝푝

푅푅

푑푑

푇푇

(4)

Table

summarizes the thermodynamic constants and other parameters used in the present study.

Thermodynamic consistency of the anelastic system requires that thermodynamic quantities are evaluated with

the reference pressure

(

) (Pauluis,

2008

). Therefore, the liquid-ice potential temperature we use is

휃휃

푙푙

(

푇푇푇푇푇

푙푙

푇푇푇

푖푖

;

푝푝

)

푇푇

(

1−

퐿퐿

푣푣푇

푇푇

푙푙

퐿퐿

푠푠푇

푇푇

푖푖

푐푐

푝푝푝푝

푇푇

)

푇

(5)

where

Π=

(

푝푝

(

푧푧

)

푝푝

)

푅푅

푚푚

∕

푐푐

푝푝푚푚

(6)

is the Exner function, evaluated with the altitude-dependent reference pressure

(

) and the constant pressure

. We take water vapor and suspended cloud condensate into account in the gas “constant“

= (1 −

)

+ (

−

)

(which is not constant because it depends on the total specific humidity

and condensate specific

humidity

) and in the isobaric specific heat

= (1 −

)

+ (

−

)

With these definitions, the anelastic governing equations in conservation form are

∇

(

휌휌

퐮퐮

)

(7)

휕휕

(

휌휌

′

퐮퐮

)

휕휕휕휕

‖Δ

⋅

(

휌휌

퐮퐮

⊗

퐮퐮

)

=−

휌휌

∇

(

훼훼

푝푝

′

)

휌휌

푏푏

퐤퐤

−

푓푓

퐤퐤

휌휌

(

퐮퐮

−

퐮퐮

푔푔

)

+∇

⋅

(

휌휌

휎휎

)

(8)

휕휕

(

휌휌

′

휃휃

푙푙

)

휕휕휕휕

+∇

⋅

(

휌휌

퐮퐮

휃휃

푙푙

)

=−

푐푐

푝푝푝푝

∇

⋅

(

휌휌

퐅퐅

푅푅

)

휌휌

푤푤

sub

휕휕휃휃

푙푙

휕휕휕휕

∇

⋅

(

휌휌

휈휈

휕휕

∇

휃휃

푙푙

)

(9)

Symbol

Name

Value

Constant reference pressure

1,000 hPa

Reference potential temperature

290 K

Gas constant of dry air

287 J (kg K)

−1

Gas constant of water vapor

461.89 J (kg K)

−1

Isobaric specific heat capacity of dry air

1004.5 J (kg K)

−1

Isobaric specific heat capacity of water vapor

1859.5 J (kg K)

−1

Specific heat capacity of liquid water

4181 J (kg K)

−1

Specific heat capacity of ice

2100 J (kg K)

−1

Specific latent heat of vapourization

2.47 MJ kg

−1

Specific latent heat of sublimation

2.83 MJ kg

−1

Freezing point temperature

273.15 K

Coriolis parameter

7.62 × 10

−5

−1

Gravitational acceleration

9.81 m s

−2

Smagorinsky constant

0.18

Turbulent Prandtl number

0.4

퐴퐴

푞푞

푡푡

Turbulent Schmidt number of water

0.4

Table 1

Thermodynamic Constants and Other Parameters Used in This Study

Journal of Advances in Modeling Earth Systems

CHAMMAS ET AL.

10.1029/2023MS003619

4 of 18

휕휕

(

휌휌

′

푞푞

푡푡

)

휕휕푡푡

+∇

⋅

(

휌휌

퐮퐮

푞푞

푡푡

)

휌휌

푤푤

sub

휕휕푞푞

푡푡

휕휕휕휕

푞푞

푡푡

∇

⋅

(

휌휌

휈휈

푡푡

∇

푞푞

푡푡

)

(10)

Here,

푏푏

푔푔

훼훼

(

휃휃

푙푙

,푞푞

푡푡

,푝푝

)

−

훼훼

(

푧푧

)

훼훼

(

푧푧

)

(11)

is the buoyancy, and

= 1/

and

= 1/

are specific volumes. The specific

volume

(

) is calculated from the approximate equation of state,

again with the reference pressure

in place of the total pressure,

훼훼

푅푅

푚푚

푇푇

푝푝

Neglected in these equations is differential settling of condensate relative to

the surrounding air and all precipitation processes.

Table

lists the variables

we use. The perturbation pressure

′ is obtained as solution to a Poisson equa

tion, which follows by taking the divergence of the momentum equation. The

numerical algorithm for solving Equations

–

is discussed in Section

2.4

2.2.

Saturation Adjustment

The temperature

and the partitioning of total water mass into the liquid

phase (specific humidity

) and ice phase (specific humidity

) are obtained

from

, and the reference pressure

by a saturation adjustment procedure

(Tao et al.,

1989

). This amounts to solving

휃휃

∗

푙푙

′

휃휃

푙푙

(12)

where

퐴퐴퐴퐴

푙푙

′

푇푇

;

푝푝

)

퐴퐴

푙푙

(

푇푇푇푇푇

∗

푙푙

푇푇푇

∗

푖푖

;

푝푝

)

is the liquid-ice potential temperature at

saturation, that is, with

푞푞

푙푙

′ −‖Δ

[

≤

,푞푞

푡푡

−

푞푞

∗

푣푣

(

푇푇,

푇푇

)

]



(

푇푇

−

푇푇

푓푓

)

(13)

and

푞푞

∗

푖푖

′ −‖Δ

[

≤

,푞푞

푡푡

−

푞푞

∗

푣푣

(

푇푇,

푇푇

)

]



(

푇푇

푓푓

−

푇푇

)

(14)

Here,

퐴퐴퐴퐴

∗

푣푣

is the saturation specific humidity, calculated as in Sridhar et al. (

2022

퐴퐴



is the Heaviside step function,

and

is the freezing point temperature. We solve the resulting nonlinear problem Equation

with the secant

method. In the presence of mixed-phase clouds, the requirement of instantaneous thermodynamic equilibrium

should be relaxed, for example, by replacing the Heaviside function in Equations

and

by a continuous phase

partitioning function (e.g., Pressel et al.,

2015

; Tao et al.,

1989

), or by carrying separate prognostic variables for

liquid and ice specific humidities. However, in the examples here we focus on warm clouds with only liquid.

2.3.

Subgrid-Scale Models

We model subgrid-scale fluxes with the turbulent viscosity model of Lilly (

1962

) and Smagorinsky (

1963

). In

this model, the turbulent viscosity is represented as

휈휈

푡푡

(

푐푐

푠푠

)

푓푓

퐵퐵

푆푆푆

(15)

where

= ‖

‖

is the 2-norm of the strain rate tensor

퐴퐴

퐒퐒

[

∇

퐮퐮

+(∇

퐮퐮

)

푇푇

]

for the resolved velocities

;

is the

Smagorinsky constant (

Table

); and Δ = (Δ

)

1/3

is the geometric mean of the grid spacings in the three

space directions. The buoyancy factor 0 ≤

≤ 1 limits the mixing length in the vertical in the case of stable

stratification; it is computed from the moist buoyancy frequency (Durran & Klemp,

1982

) as described in Pressel

et al. (

2017

). The diffusivities of the liquid-ice potential temperature and total specific humidity are obtained

from the turbulent viscosity

by division by constant turbulent Prandtl and Schmidt numbers (Table

Variable

Definition

Units

Density of moist air

kg m

−3

Specific volume of moist air

−1

Velocity of moist air

m s

−1

Prescribed geostrophic velocity

m s

−1

sub

Prescribed subsidence velocity

m s

−1

Pressure

Buoyancy

m s

−2

Vertical unit vector

Temperature

Specific gas “constant” of moist air

J kg

−1

Isobaric specific heat of moist air

J kg

−1

Subgrid-scale stress per unit mass

−2

Radiative energy flux

W m kg

−1

Total water specific humidity

kg/kg

Water vapor specific humidity

kg/kg

Liquid water specific humidity

kg/kg

Ice specific humidity

kg/kg

Turbulent viscosity

−1

Altitude

Table 2

Definitions of Variables

Journal of Advances in Modeling Earth Systems

CHAMMAS ET AL.

10.1029/2023MS003619

5 of 18

To emulate a radiation condition at the upper boundary, we include a sponge layer that occupies the top 5% of

the domain and absorbs upward propagating waves. The sponge is implemented as a linear Rayleigh damping

layer (Durran & Klemp,

1983

), in which the horizontal velocity is relaxed toward the geostrophic wind velocity

and the vertical velocity is relaxed to zero. To avoid reflections at the interface between the sponge layer and the

undamped flow outside, we use a relaxation coefficient that ensures a gradual onset of the sponge layer (Klemp

& Lilly,

1978

), reaching 0.25 s

−1

at the top of the domain.

2.4.

Numerical Solution

We discretize the governing equations with the finite-difference method. All discrete operators are expressed on

a collocated mesh. All diffusion terms are computed with a second-order central difference scheme. The advec

tion terms in Equations

–

are discretized with the third-order QUICK (Quadratic Upstream Interpolation for

Convective Kinematics) scheme. While the QUICK scheme is upwind-biased, it is not monotonicity preserving.

This implies that in interaction with subgrid-scale diffusion, it can create spurious mixing in regions of sharp

gradients, for example, at a sharp inversion topping a boundary layer, with potentially deleterious effects on the

simulation of stratocumulus clouds (Bretherton et al.,

1999

; Pressel et al.,

2017

). As we will see, these effects are

real but minor in our case. Alternatively, one may construct a monotone version of the QUICK scheme by apply

ing flux limiters (Zalesak,

1979

; D. E. Stevens & Bretherton,

1996

). To test the effects of the advection scheme

near the inversion, we also implemented a third-order WENO scheme, reconstructing the advective fluxes on cell

faces and using a Lax-Friedrichs numerical flux, as in standard finite-volume methods.

An explicit iterative scheme (Wang et al.,

2022

) is employed for the time advancement of the numerical solutions.

This scheme provides an iterative representation to the Crank-Nicolson method, which avoids the computational

complexity of solving a high-dimensional linear system of equations. Specifically, the momentum Equation

solved with a predictor-corrector approach. At the prediction step of sub-iteration

+ 1, the momentum equation

is solved in discrete form as

휌휌

퐮퐮

−

(

휌휌

퐮퐮

)

푛푛

푡푡

=−

휌휌

∇

(

훼훼

푝푝

′

푘푘

)

퐑퐑

푛푛

(16)

with

퐑퐑

=−∇

⋅

[

(

휌휌

퐮퐮

)

푛푛

⊗

퐮퐮

푛푛

]

+∇

⋅

(

휌휌

휎휎

푛푛

)

휌휌

푏푏

푛푛

퐤퐤

−

푓푓

퐤퐤

휌휌

(

퐮퐮

푛푛

−

퐮퐮

푔푔

)

(17)

where

퐴퐴

(

⋅

)

denotes a prediction of a variable at step

+ 1 in sub-iteration

+ 1; (⋅)

is the solution of a variable at

step

+ 1 obtained from sub-iteration

. Variables at state

퐴퐴

(

⋅

)

푛푛

are estimated as

퐴퐴

(

⋅

)

푛푛

[

(

⋅

)

푘푘

⋅

)

푛푛

]

∕2

. Note

that the prediction of the momentum

퐴퐴

휌휌

퐮퐮

in sub-iteration

+ 1 is evaluated with the pressure from the previous

sub-iteration. The correct momentum

퐴퐴

(

휌휌

퐮퐮

)

푘푘

needs to be computed with the pressure at sub-iteration

+ 1,

which can be expressed similarly to Equation

(

퐮퐮

)

푘푘

−

(

휌휌

퐮퐮

)

푛푛

푡푡

=−

휌휌

∇

(

훼훼

푝푝

′

푘푘

)

퐑퐑

푛푛

(18)

Subtracting Equation

from Equation

yields

(

퐮퐮

)

푘푘

−

휌휌

퐮퐮

푡푡

=−

휌휌

∇

(

훼훼

푝푝

′

푘푘

−

훼훼

푝푝

′

푘푘

)

=−

휌휌

∇

(

훼훼

훿훿푝푝

)

(19)

where

δp

′

−

′

is the pressure correction from sub-iteration

+ 1.

Taking the divergence of Equation

and applying mass conservation at sub-iteration

+ 1 leads to a generalized

Poisson equation for the pressure correction:

∇

(

훼훼

훿훿훿훿

)

훼훼

푡푡

[

∇

⋅

(

휌휌

퐮퐮

)

−∇

⋅

(

휌휌

퐮퐮

)

푘푘

]

훼훼

푡푡

∇

⋅

(

휌휌

퐮퐮

)

(20)

To ensure numerical consistency and eliminate the checkerboard effect due to the collocated mesh representation,

we introduce an additional correction term when solving Equation

, which is described in Appendix