Evolution of a Buoyancy Gradient in 2D Boussinesq Framework 1 ...
Evolution of a Buoyancy Gradient in 2D Boussinesq Framework Benjamin Mayersohn Courant Institute of Mathematical Sciences
1
Introduction
The Boussinesq equations are a simplification to the full 3D Euler equations for a fluid with small density variations, such as a stratified ocean. Consider the original Euler equations in a rotating, stratified fluid, with no heat source (here cs is the speed of sound): ρ
D~v + ρf kˆ × ~v = −∇p − ρg kˆ Dt Dρ + ρ∇ · ~v = 0 Dt 1 Dp Dρ − 2 =0 Dt cs Dt
Assuming small density variations, we can proceed like Vallis [1]: let ρ = ρ0 + δρ(x, y, z, t) = ρ0 + ρˆ(z) + ρ0 (x, y, z, t), and assume that the density variations are insignificant except when they are scaled by gravity terms, which are relatively large. We also define a reference pressure that relates to the large-scale density 0 in the following way: ∂p ∂z = −ρ0 g, which is a standard assumption when hydrostatic balance holds (or when the flow is steady in the z-direction). Substituting these expressions into the above equations and making a few more approximations based on the relative sizes of the terms, we arrive at the equations (getting rid of stylistic ticks)
D~u 1 + f kˆ × ~u = − ∇h p Dt ρ0 Dw 1 ∂p =− +b Dt ρ0 ∂z Db + N 2w = 0 Dt ∇ · ~v = 0 0
2 where ∇h is the horizontal gradient, b0 = − gρ ρ0 and N =
1.1
∂ˆ b ∂z
=
∂ ∂z
− gρρ0ˆ .
2D Flows
The Boussinesq equations allow for flows that vary in the vertical and only one horizontal dimension, given that the initial condition is also independent of that other dimension. Say we want a solution without
1
y-dependence. Then ∂u ∂t ∂v ∂t ∂w ∂t ∂b ∂t ux + wz
= −∂x φ + f v − uux − wuz = −f u − uvx − wvz = −∂z φ − uwx − wwz + b = −ubx − wbz − N 2 w =0
where φ = ρp0 . Such flows admit a streamfunction ψ, defined as the Laplacian of the y-component of vorticity, ω y = uz − wx . We can reduce the system to three dynamical equations for buoyancy, the meridional velocity, and the streamfunction (plus the elliptic equation that relates vorticity and the streamfunction). The streamfunction is defined as u = ∂z ψ w = −∂x ψ so that uz − wx = ω y = ∇2 ψ. Incompressibility of course holds. The dynamical equations are then ∂∇2 ψ + J(∇2 ψ, ψ) = f vz − bx ∂t ∂b + J(b, ψ) = N 2 ψx ∂t ∂v + J(v, ψ) = −f ψz ∂t where J(A, B) = Ax Bz − Az Bx is the Jacobian in x-z space.
2
Buoyancy-Induced Flow in 2D
Though flows in 2D are relatively simple, the equations are still interesting enough to consider a variety of situations. One in which I am interested in is the structure and evolution of surface fronts, which are ubiquitous throughout the ocean [2]. To study the evolution of a temperature gradient one may start with the buoyancy equation: Db + N 2 w = bt + ubx + wbz + N 2 w = 0 Dt Say our anomalous buoyancy gradient changes mainly in the x-direction. Then we can take the x-derivative of both sides: (bx )t + (~u)x · ∇b + ~u · ∇bx + N 2 wx = 0 Dbx = −(~u)x · ∇b − N 2 wx Dt The nonlinear term −(~u)x · ∇b often dominates and enhances gradients if they are aligned favorably with a strain field (e.g., if a strain field converges on an axis parallel to the gradient). Say the buoyancy anomaly 2
gradient is zero in the z-direction (note that there is still a total buoyancy gradient, but it is the deviation from hydrostatic balance that we are considering) and we have a flow that converges in x (ux < 0), with no flow in the vertical direction. Then Dbx = −ux bx Dt which indicates that the magnitude bx would grow exponentially, at least early on. In addition, consider the momentum equations, and write the equations for wx and uz : Dwx = −~ux · ∇w − ∂z (φx ) + bx Dt Duz = −~uz · ∇u + f vz − ∂z (φx ) Dt Dω y D =⇒ = (uz − wx ) = ~ux · ∇w − ~uz · ∇(u) + f vz − bx Dt Dt = ux wx + wx wz − uz ux − wz uz + f vz − bx (((( (((( (−uz( (( (−w w( =( z ) − wz uz ) + f vz − bx = f vz − bx x + wx wz ) + ( (z( ( (−w If the flow were in geostrophic balance (so that f vz = bx , which is often the case in large-scale geophysical regimes where rotational effects are dominant), then the vorticity would be materially conserved. Say we initially have no flow and a negative buoyancy anomaly (i.e. positive density anomaly) in the center of a channel (x-periodic, rigid walls in c) near the surface, which decays to zero at the edges of the periodic domain (see Figure 1). Then we begin out of balance, and the vorticity must adjust. bx > 0 to the right of the anomaly and bx < 0 at the left. From the y-vorticity evolution equation we expect a clockwise rotational response to the left of the anomaly and a counter-clockwise response to the right. This double-celled gyre has the effect of subducting the negatively-buoyant anomaly and replacing it with warmer water at the surface, thus restratifying the water column.
Figure 1: Initial surface buoyancy gradient. Distribution exponentially decays in z to the bottom The plots in Figure 2 show the evolution of the velocity field and buoyancy contours in a simulation run using the model I am about to describe, and confirm the conclusion we have come to mathematically. 3
(a)
(b)
(c)
(d)
Figure 2: Buoyancy contour ((a) and (c)) and velocity field ((b) and (d)) plots for initial surface buoyancy gradient described in Figure 1. Taken at times t=10 and t=40, for a simulation that runs for 100 seconds total.
3
The Projection Method
We will not take advantage of the streamfunction formulation posed earlier. We will instead use a classic method by Chorin [3], and follow the implementation considered by Durran [4]. The projection method begins by integrating the momentum equations in time: n+1 tZ
∂v dt = − ∂t
tn
n+1 tZ
n+1 tZ
F~~v (~v n , bn )dt
∇φdt + tn
tn
˜ The next step is to define an “average” pressure gradient ∇φ: ˜n+1
∇φ
n+1 tZ
∇φdt
∆t = tn
4
Substituting this back into the equations, and rewriting the integral of ~vt as the difference between the velocities at the two times, we get ~v n+1 − ~v n = −∆t∇φ˜n+1 +
n+1 tZ
F~~v (~v n , bn )dt
tn
Now we can define an intermediate velocity v˜ and substitute it back into the previous equation ~v˜ = ~v n +
n+1 tZ
F~~v (~v n , bn )dt
tn
~v
n+1
= ~v˜ − ∆t∇φ˜n+1
These two equations are all we need to construct the pressure method. The first equation may be solved numerically via some time stepping scheme. The second equation becomes an elliptic equation for pressure if we take the divergence of it. This is because ∇ · ~v n+1 = 0, since we must enforce continuity of the velocity field at each time step. Thus ∇2 φ˜n+1 =
∇ · ~v˜ ∆t
We can also derive the boundary condition for a rigid boundary from that same equation. If we take the dot product between each side of the equation and the normal to the boundary, and we enforce the notion that the component of the true velocity normal to the boundary must be zero, then ∇φ˜n+1 · n ˆ=
~v˜ · n ˆ ∆t
Thus the three phases involved are 1. Solve the equation for the intermediate velocity: ~v˜ = ~v n +
tn+1 R
F~~v (~v n , bn )dt
tn
2. Solve for the pressure and use the appropriate boundary conditions: ∇2 φ˜n+1 =
˜ ∇·~ v ∆t
3. Update the velocity: ~v n+1 = ~v˜ − ∆t∇φ˜n+1
4
Methods + Implementation
For my project, I implemented the situation described in Figures 1-2 (i.e. the free evolution of a buoyancy gradient with no initial flow). The initial condition is given by the equation: √ z−Lz f π x − Lx /2 + d x − Lx /2 − d b= us erf − erf e h0 2 l0 l0 where l0 , h0 , d are parameters controlling the horizontal and vertical decay rates of the buoyancy anomaly, Lx , Lz are the x and z extents of the domain, us is a velocity scale (used to make the units of buoyancy, which are units of acceleration, consistent) and f is the Coriolis parameter, which is approximately 1 × 10−4 /s.
5
The domain was constructed as a 500 m by 100 m box, with periodic boundary conditions in x and free-slip conditions in z (uz = vz = 0). The x boundary conditions are of little consequence, as the box is wide enough such that boundary contributions are minimal. nz , where nz is the number of grid points in the z-direction The domain is divided into z-stacks of size np and np is the number of processors. The domain is not parallelized in x. The problem is parallelized using MPICH2, which utilizes MPI (Message-Passing Interface) to distribute memory in parallel between multiple cores. The model was run on New York University’s Crunchy 3 CIMS cluster, which consists of four 8-core AMD Opteron 6136 2.4GHz processors.
Spatial derivatives are calculated using second order finite differences, in both x and z. The periodic nature of x makes Fourier methods better suited to solve the problem in this direction, but for the sake of time I opted to use finite differences. Phase 1 (solving for the intermediate velocity) is implemented using a leap-frog time stepper with a RobertAsselin filter [4], after an initial Euler step and a secondary unfiltered leap frog step. The Robert coefficient is assigned to 0.06, as is common for problems in atmosphere-ocean sciences. Phase 2 (solving for the pressure) is implemented using a Jacobi solver. This is quite slow, converging like O(N 2 ), without the help of a method such as multigrid that converges like O(N ). This is likely the step that makes reduces the effectiveness of parallelization. Phase 3 (updating the velocity) requires no special methods, and is simply called by combining the terms calculated in the previous two steps.
4.1
Steps, Methods + MPI Calls
Listed below are the steps taken for each time step: 1. MPI pass: exchange u,v,w,φ,b between processors 2. Calculate derivatives found on RHS of evolution equations for u,v,w,b 3. Pass RHS to time stepper, solve for intermediate quantities 4. MPI pass: Exchange w, ˜ which is needed for ∇ · ~v˜ 5. Calculate divergence (MPI pass for z-derivatives), set as RHS for elliptic equation for pressure 6. Solve for interior values via Jacobi (uses MPI again for z-derivatives) + update boundary values for pressure using BC’s 7. Calculate pressure gradients (uses MPI for z-derivatives) and use these to calculate velocity at t+dt MPI is called quite often using this strategy, which will inevitably create overhead, particularly in the Jacobi step. All exchanges are made between neighboring processors, which is all that is required for finite difference derivatives in the z-direction, using MPI ISend and MPI IRecv, followed by an MPI Wait to confirm that the nonblocking requests have been completed by each processor.
5
Results
The model was run for a variety of different domain lengths and processors, with the hopes that parallelization would lead to a speed increase. The results are included in Table 1. In the first box, the model is run with an equal number of points in x and z. We see an inconsistent trend, with a speedup between 2 and 4 processors,
6
followed by a subsequent slowdown. This implies that the point at which overhead overtakes the value of parallelization is somewhere between 4 and 8 processors. In the second box, we increase the number of points in z and decrease the number of points in x. We see here that any kind of parallelization actually slows down the computations, as increasing the processors increases the number of exchanges in the z-direction, which there are more of because of the increase in the number of z-points. The third box demonstrates that overhead is truly the issue at hand. We see that when we have more points in x than in z, adding more processors results in a computational speedup. The reason is that xdifferentiation is an embarrassingly parallel operation, as there is no dependence between processors. Adding more processors simply reduces the number of x-derivatives per processor that need to be computed. More strain is placed on z-derivatives, but there are fewer points in z here, so the overall increase in effort here is undermined by the relative ease of taking derivatives in the x-direction. Fortunately problems in atmosphereocean science often require greater resolution in x than in z, since the aspect ratio, or the ratio between the z and x extents of the domain, is often small. nz=128, nx=128, dt=0.1, T=100 np Time (s) 2 16.54 4 11.32 8 15.36 16 19.32
nz=128, nx=32, dt=0.1, T=100 np Time (s) 2 3.56 4 4.49 8 7.87
nz=128, nx=256, dt=0.1, T=10 np Time (s) 2 3.4 4 2.34 8 2.07 16 1.8 Table 1: Runs for various domain sizes/processors. Here, nz and nx are the number of z and x grid points (respectively), dt is the time step and T is the end time.
6
Summary + Conclusion
We have implemented the projection method for a the 2D Boussinesq equations, with an initial negative buoyancy gradient at the surface of the domain. We parallelized in the z-direction by using MPI to distribute the memory across processors and solved for a variety of gridpoint distributions. We found that the implementation is quite inefficient, other than in the case where the most expensive operations can be made embarrassingly parallel (i.e. taking x-derivatives). The cause is likely the use of a Jacobi solver, which has smooth but slow convergence properties. Using a multigrid method would likely allow for speedup, and I hope to implement this in the future when testing more complicated flows.
References [1] Vallis, G.K. (2006). Atmospheric and Oceanic Fluid Dynamics. 1st Edition: Cambridge University Press [2] McWilliams, J.C., Colas, F., & Molemaker, M. J. (2009). Cold filamentary intensification and oceanic surface convergence lines. Geophysical Research Letters, 36. 7
[3] Chorin, A.J. (1967). The numerical solution of the Navier-Stokes equations for an incompressible fluid. Bull. Am. Math. Soc. 73: 928-931 [4] Durran, D.R. (2010). Numerical Methods for Fluid Dynamics: With Applications to Geophysics (Texts in Applied Mathematics. 2nd Edition: Springer.
8
1
Introduction
The Boussinesq equations are a simplification to the full 3D Euler equations for a fluid with small density variations, such as a stratified ocean. Consider the original Euler equations in a rotating, stratified fluid, with no heat source (here cs is the speed of sound): ρ
D~v + ρf kˆ × ~v = −∇p − ρg kˆ Dt Dρ + ρ∇ · ~v = 0 Dt 1 Dp Dρ − 2 =0 Dt cs Dt
Assuming small density variations, we can proceed like Vallis [1]: let ρ = ρ0 + δρ(x, y, z, t) = ρ0 + ρˆ(z) + ρ0 (x, y, z, t), and assume that the density variations are insignificant except when they are scaled by gravity terms, which are relatively large. We also define a reference pressure that relates to the large-scale density 0 in the following way: ∂p ∂z = −ρ0 g, which is a standard assumption when hydrostatic balance holds (or when the flow is steady in the z-direction). Substituting these expressions into the above equations and making a few more approximations based on the relative sizes of the terms, we arrive at the equations (getting rid of stylistic ticks)
D~u 1 + f kˆ × ~u = − ∇h p Dt ρ0 Dw 1 ∂p =− +b Dt ρ0 ∂z Db + N 2w = 0 Dt ∇ · ~v = 0 0
2 where ∇h is the horizontal gradient, b0 = − gρ ρ0 and N =
1.1
∂ˆ b ∂z
=
∂ ∂z
− gρρ0ˆ .
2D Flows
The Boussinesq equations allow for flows that vary in the vertical and only one horizontal dimension, given that the initial condition is also independent of that other dimension. Say we want a solution without
1
y-dependence. Then ∂u ∂t ∂v ∂t ∂w ∂t ∂b ∂t ux + wz
= −∂x φ + f v − uux − wuz = −f u − uvx − wvz = −∂z φ − uwx − wwz + b = −ubx − wbz − N 2 w =0
where φ = ρp0 . Such flows admit a streamfunction ψ, defined as the Laplacian of the y-component of vorticity, ω y = uz − wx . We can reduce the system to three dynamical equations for buoyancy, the meridional velocity, and the streamfunction (plus the elliptic equation that relates vorticity and the streamfunction). The streamfunction is defined as u = ∂z ψ w = −∂x ψ so that uz − wx = ω y = ∇2 ψ. Incompressibility of course holds. The dynamical equations are then ∂∇2 ψ + J(∇2 ψ, ψ) = f vz − bx ∂t ∂b + J(b, ψ) = N 2 ψx ∂t ∂v + J(v, ψ) = −f ψz ∂t where J(A, B) = Ax Bz − Az Bx is the Jacobian in x-z space.
2
Buoyancy-Induced Flow in 2D
Though flows in 2D are relatively simple, the equations are still interesting enough to consider a variety of situations. One in which I am interested in is the structure and evolution of surface fronts, which are ubiquitous throughout the ocean [2]. To study the evolution of a temperature gradient one may start with the buoyancy equation: Db + N 2 w = bt + ubx + wbz + N 2 w = 0 Dt Say our anomalous buoyancy gradient changes mainly in the x-direction. Then we can take the x-derivative of both sides: (bx )t + (~u)x · ∇b + ~u · ∇bx + N 2 wx = 0 Dbx = −(~u)x · ∇b − N 2 wx Dt The nonlinear term −(~u)x · ∇b often dominates and enhances gradients if they are aligned favorably with a strain field (e.g., if a strain field converges on an axis parallel to the gradient). Say the buoyancy anomaly 2
gradient is zero in the z-direction (note that there is still a total buoyancy gradient, but it is the deviation from hydrostatic balance that we are considering) and we have a flow that converges in x (ux < 0), with no flow in the vertical direction. Then Dbx = −ux bx Dt which indicates that the magnitude bx would grow exponentially, at least early on. In addition, consider the momentum equations, and write the equations for wx and uz : Dwx = −~ux · ∇w − ∂z (φx ) + bx Dt Duz = −~uz · ∇u + f vz − ∂z (φx ) Dt Dω y D =⇒ = (uz − wx ) = ~ux · ∇w − ~uz · ∇(u) + f vz − bx Dt Dt = ux wx + wx wz − uz ux − wz uz + f vz − bx (((( (((( (−uz( (( (−w w( =( z ) − wz uz ) + f vz − bx = f vz − bx x + wx wz ) + ( (z( ( (−w If the flow were in geostrophic balance (so that f vz = bx , which is often the case in large-scale geophysical regimes where rotational effects are dominant), then the vorticity would be materially conserved. Say we initially have no flow and a negative buoyancy anomaly (i.e. positive density anomaly) in the center of a channel (x-periodic, rigid walls in c) near the surface, which decays to zero at the edges of the periodic domain (see Figure 1). Then we begin out of balance, and the vorticity must adjust. bx > 0 to the right of the anomaly and bx < 0 at the left. From the y-vorticity evolution equation we expect a clockwise rotational response to the left of the anomaly and a counter-clockwise response to the right. This double-celled gyre has the effect of subducting the negatively-buoyant anomaly and replacing it with warmer water at the surface, thus restratifying the water column.
Figure 1: Initial surface buoyancy gradient. Distribution exponentially decays in z to the bottom The plots in Figure 2 show the evolution of the velocity field and buoyancy contours in a simulation run using the model I am about to describe, and confirm the conclusion we have come to mathematically. 3
(a)
(b)
(c)
(d)
Figure 2: Buoyancy contour ((a) and (c)) and velocity field ((b) and (d)) plots for initial surface buoyancy gradient described in Figure 1. Taken at times t=10 and t=40, for a simulation that runs for 100 seconds total.
3
The Projection Method
We will not take advantage of the streamfunction formulation posed earlier. We will instead use a classic method by Chorin [3], and follow the implementation considered by Durran [4]. The projection method begins by integrating the momentum equations in time: n+1 tZ
∂v dt = − ∂t
tn
n+1 tZ
n+1 tZ
F~~v (~v n , bn )dt
∇φdt + tn
tn
˜ The next step is to define an “average” pressure gradient ∇φ: ˜n+1
∇φ
n+1 tZ
∇φdt
∆t = tn
4
Substituting this back into the equations, and rewriting the integral of ~vt as the difference between the velocities at the two times, we get ~v n+1 − ~v n = −∆t∇φ˜n+1 +
n+1 tZ
F~~v (~v n , bn )dt
tn
Now we can define an intermediate velocity v˜ and substitute it back into the previous equation ~v˜ = ~v n +
n+1 tZ
F~~v (~v n , bn )dt
tn
~v
n+1
= ~v˜ − ∆t∇φ˜n+1
These two equations are all we need to construct the pressure method. The first equation may be solved numerically via some time stepping scheme. The second equation becomes an elliptic equation for pressure if we take the divergence of it. This is because ∇ · ~v n+1 = 0, since we must enforce continuity of the velocity field at each time step. Thus ∇2 φ˜n+1 =
∇ · ~v˜ ∆t
We can also derive the boundary condition for a rigid boundary from that same equation. If we take the dot product between each side of the equation and the normal to the boundary, and we enforce the notion that the component of the true velocity normal to the boundary must be zero, then ∇φ˜n+1 · n ˆ=
~v˜ · n ˆ ∆t
Thus the three phases involved are 1. Solve the equation for the intermediate velocity: ~v˜ = ~v n +
tn+1 R
F~~v (~v n , bn )dt
tn
2. Solve for the pressure and use the appropriate boundary conditions: ∇2 φ˜n+1 =
˜ ∇·~ v ∆t
3. Update the velocity: ~v n+1 = ~v˜ − ∆t∇φ˜n+1
4
Methods + Implementation
For my project, I implemented the situation described in Figures 1-2 (i.e. the free evolution of a buoyancy gradient with no initial flow). The initial condition is given by the equation: √ z−Lz f π x − Lx /2 + d x − Lx /2 − d b= us erf − erf e h0 2 l0 l0 where l0 , h0 , d are parameters controlling the horizontal and vertical decay rates of the buoyancy anomaly, Lx , Lz are the x and z extents of the domain, us is a velocity scale (used to make the units of buoyancy, which are units of acceleration, consistent) and f is the Coriolis parameter, which is approximately 1 × 10−4 /s.
5
The domain was constructed as a 500 m by 100 m box, with periodic boundary conditions in x and free-slip conditions in z (uz = vz = 0). The x boundary conditions are of little consequence, as the box is wide enough such that boundary contributions are minimal. nz , where nz is the number of grid points in the z-direction The domain is divided into z-stacks of size np and np is the number of processors. The domain is not parallelized in x. The problem is parallelized using MPICH2, which utilizes MPI (Message-Passing Interface) to distribute memory in parallel between multiple cores. The model was run on New York University’s Crunchy 3 CIMS cluster, which consists of four 8-core AMD Opteron 6136 2.4GHz processors.
Spatial derivatives are calculated using second order finite differences, in both x and z. The periodic nature of x makes Fourier methods better suited to solve the problem in this direction, but for the sake of time I opted to use finite differences. Phase 1 (solving for the intermediate velocity) is implemented using a leap-frog time stepper with a RobertAsselin filter [4], after an initial Euler step and a secondary unfiltered leap frog step. The Robert coefficient is assigned to 0.06, as is common for problems in atmosphere-ocean sciences. Phase 2 (solving for the pressure) is implemented using a Jacobi solver. This is quite slow, converging like O(N 2 ), without the help of a method such as multigrid that converges like O(N ). This is likely the step that makes reduces the effectiveness of parallelization. Phase 3 (updating the velocity) requires no special methods, and is simply called by combining the terms calculated in the previous two steps.
4.1
Steps, Methods + MPI Calls
Listed below are the steps taken for each time step: 1. MPI pass: exchange u,v,w,φ,b between processors 2. Calculate derivatives found on RHS of evolution equations for u,v,w,b 3. Pass RHS to time stepper, solve for intermediate quantities 4. MPI pass: Exchange w, ˜ which is needed for ∇ · ~v˜ 5. Calculate divergence (MPI pass for z-derivatives), set as RHS for elliptic equation for pressure 6. Solve for interior values via Jacobi (uses MPI again for z-derivatives) + update boundary values for pressure using BC’s 7. Calculate pressure gradients (uses MPI for z-derivatives) and use these to calculate velocity at t+dt MPI is called quite often using this strategy, which will inevitably create overhead, particularly in the Jacobi step. All exchanges are made between neighboring processors, which is all that is required for finite difference derivatives in the z-direction, using MPI ISend and MPI IRecv, followed by an MPI Wait to confirm that the nonblocking requests have been completed by each processor.
5
Results
The model was run for a variety of different domain lengths and processors, with the hopes that parallelization would lead to a speed increase. The results are included in Table 1. In the first box, the model is run with an equal number of points in x and z. We see an inconsistent trend, with a speedup between 2 and 4 processors,
6
followed by a subsequent slowdown. This implies that the point at which overhead overtakes the value of parallelization is somewhere between 4 and 8 processors. In the second box, we increase the number of points in z and decrease the number of points in x. We see here that any kind of parallelization actually slows down the computations, as increasing the processors increases the number of exchanges in the z-direction, which there are more of because of the increase in the number of z-points. The third box demonstrates that overhead is truly the issue at hand. We see that when we have more points in x than in z, adding more processors results in a computational speedup. The reason is that xdifferentiation is an embarrassingly parallel operation, as there is no dependence between processors. Adding more processors simply reduces the number of x-derivatives per processor that need to be computed. More strain is placed on z-derivatives, but there are fewer points in z here, so the overall increase in effort here is undermined by the relative ease of taking derivatives in the x-direction. Fortunately problems in atmosphereocean science often require greater resolution in x than in z, since the aspect ratio, or the ratio between the z and x extents of the domain, is often small. nz=128, nx=128, dt=0.1, T=100 np Time (s) 2 16.54 4 11.32 8 15.36 16 19.32
nz=128, nx=32, dt=0.1, T=100 np Time (s) 2 3.56 4 4.49 8 7.87
nz=128, nx=256, dt=0.1, T=10 np Time (s) 2 3.4 4 2.34 8 2.07 16 1.8 Table 1: Runs for various domain sizes/processors. Here, nz and nx are the number of z and x grid points (respectively), dt is the time step and T is the end time.
6
Summary + Conclusion
We have implemented the projection method for a the 2D Boussinesq equations, with an initial negative buoyancy gradient at the surface of the domain. We parallelized in the z-direction by using MPI to distribute the memory across processors and solved for a variety of gridpoint distributions. We found that the implementation is quite inefficient, other than in the case where the most expensive operations can be made embarrassingly parallel (i.e. taking x-derivatives). The cause is likely the use of a Jacobi solver, which has smooth but slow convergence properties. Using a multigrid method would likely allow for speedup, and I hope to implement this in the future when testing more complicated flows.
References [1] Vallis, G.K. (2006). Atmospheric and Oceanic Fluid Dynamics. 1st Edition: Cambridge University Press [2] McWilliams, J.C., Colas, F., & Molemaker, M. J. (2009). Cold filamentary intensification and oceanic surface convergence lines. Geophysical Research Letters, 36. 7
[3] Chorin, A.J. (1967). The numerical solution of the Navier-Stokes equations for an incompressible fluid. Bull. Am. Math. Soc. 73: 928-931 [4] Durran, D.R. (2010). Numerical Methods for Fluid Dynamics: With Applications to Geophysics (Texts in Applied Mathematics. 2nd Edition: Springer.
8
