Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Avoid double staggering of ww in vertical advection of geopotential #1338

Merged
merged 2 commits into from
Feb 25, 2021
Merged

Avoid double staggering of ww in vertical advection of geopotential #1338

merged 2 commits into from
Feb 25, 2021

Conversation

matzegoebel
Copy link
Contributor

@matzegoebel matzegoebel commented Dec 15, 2020 •
edited by davegill
Loading

TYPE: enhancement

KEYWORDS: geopotential,vertical advection

SOURCE: Matthias Göbel (University of Innsbruck)

DESCRIPTION OF CHANGES:
The geopotential equation in WRF is formulated in the following (advective) form where staggering/destaggering is
marked with overbars:

The vertical advection term thus employs a double staggering of Omega. First, Omega is destaggered to mass levels
and multiplied with the phi gradient. Then the product is staggered back to w-levels. This formulation is analogous to the horizontal terms, but the double averaging of Omega may introduce additional inaccuracies and is not necessary.

The phi gradient can be first staggered to w-levels and then multiplied with the unaveraged Omega:

For consistency and to avoid numerical instabilities, the same thing has to be done in the acoustic step.

The new namelist option phi_adv_z allows switching between the two formulations:
phi_adv_z = 1: old formulation (default)
phi_adv_z = 2: new formulation

LIST OF MODIFIED FILES:
Registry/Registry.EM_COMMON
dyn_em/module_big_step_utilities_em.F
dyn_em/module_small_step_em.F
run/README.namelist

TESTS CONDUCTED:
Idealized test cases em_les and em_hill2d_x.
The vertical velocity is affected by the modification. Updrafts and downdrafts in the em_hill2d_x case become ~10%
stronger with the new formulation:
W_xz
In the LES case stronger updrafts and downdrafts occur more frequently:
hist_w

I came across this issue when trying to transform advective components of temperature, humidity, and momentum into
Cartesian form in a consistent way, using w for the vertical advection. Without the proposed modifications I could not
close the budget, because the vertical advection of phi was not consistent with the vertical advection in the other budget
equations.

RELEASE NOTE: A new namelist option is included that allows a user to avoiding double stagger averaging of Omega in the vertical advection of geopotential. The new namelist option phi_adv_z allows switching between the two formulations: phi_adv_z = 1 (the old formulation, which remains the default), and phi_adv_z = 2 (the new formulation).

@davegill
Copy link
Contributor

@dudhia @weiwangncar
Folks,
Would you guys review this please

@dudhia
Copy link
Collaborator

dudhia commented Dec 29, 2020 via email

@dudhia
Copy link
Collaborator

dudhia commented Dec 29, 2020

I will say it looks OK and this discretization is not detailed in the ARW Tech Note, so nothing to modify there.

@weiwangncar
Copy link
Collaborator

@skamaroc Bill, can you take a look at this PR? Thanks!

@dudhia
Copy link
Collaborator

dudhia commented Jan 12, 2021

@matzegoebel can you indicate how much effect this has on results? Is it random noise or systematic in the diffs?

@dudhia
Copy link
Collaborator

dudhia commented Jan 12, 2021

Also would be helpful to show the change in discretization as an equation (not sure how but something like)
overbar [ overbar(w) (d phi d eta) ]
changed to
w overbar (d phi d eta)
where overbar is the vertical average.

@davegill
Copy link
Contributor

@dudhia @matzegoebel
Jimy,
Of all of the WRF contributors, you have asked the right person to show you the discretized equations!

@dudhia
Copy link
Collaborator

dudhia commented Jan 12, 2021 via email

@matzegoebel
Copy link
Contributor Author

I updated the commit message to include the corresponding equations.

@matzegoebel
Copy link
Contributor Author

matzegoebel commented Jan 13, 2021 •
edited
Loading

I ran the em_les test case at 1 min output interval with and without the modifications . Here are the histogram plots of the geopotential difference at four different heights:
hist_ph

Judging from these plots, the differences seem to be random rather than systematic. Or do you have something else in mind that I should look at?

@matzegoebel
Copy link
Contributor Author

matzegoebel commented Jan 13, 2021 •
edited
Loading

When averaging over x and y the differences become quite small:
scatter_ph
The plot shows the values for all times and heights.

@dudhia
Copy link
Collaborator

dudhia commented Jan 13, 2021

Thanks for adding the equations. I was just thinking of an ncdiff plot to show a sample effect which would illustrate noise versus systematic change. If there is any evidence that the change does not introduce noise, that would be great too, since it now uses a less vertically averaged w.

@matzegoebel
Copy link
Contributor Author

matzegoebel commented Jan 13, 2021 •
edited
Loading

I'm not quite sure what you mean exactly. Here are instantaneous time-height difference plots for different xy locations:
pcolor_ph
I cannot really see much in there. Only that there is some feature at bottom_top_stag=20 which is the boundary layer height.
And here is a kernel density estimation of the perturbation from the xy-average for data points between 30 and 60 min integration and for the lowest 20 gridpoints in height:
hist_all
The distributions are almost identical, so I think there should not be more noise (i.e. broader PDF) with the new version.

@dudhia
Copy link
Collaborator

dudhia commented Jan 13, 2021 via email

@dudhia dudhia changed the base branch from release-v4.2.2 to develop January 13, 2021 20:51
@dudhia dudhia requested review from a team as code owners January 13, 2021 20:51
@dudhia
Copy link
Collaborator

dudhia commented Jan 13, 2021

I have moved this to a develop branch PR and removed the bug label

@matzegoebel
Copy link
Contributor Author

Oh yes that's right, it's not actually a bug. I also changed the base branch of my PR.

@matzegoebel
Copy link
Contributor Author

Sure! Here is the full time-series of the absolute mu change averaged over xy and a close-up on the first 10 minutes:
dmu_timeseries
dmu_timeseries_start

@matzegoebel
Copy link
Contributor Author

it's the em_les test case

@dudhia
Copy link
Collaborator

dudhia commented Jan 20, 2021 via email

@davegill
Copy link
Contributor

davegill commented Feb 9, 2021

@dudhia @weiwangncar @matzegoebel
Folks,
Is there a status on this yet? Matthias has provided about a zillion plots and equations.

@dudhia
Copy link
Collaborator

dudhia commented Feb 9, 2021 via email

@dudhia
Copy link
Collaborator

dudhia commented Feb 11, 2021

@matzegoebel I have heard back from Bill Skamarock who wants to be cautious but he is OK with this being added via a switch, so that people can test it more for themselves before we commit to it. I think a sample difference xy plot of fields at some time into the simulation may also be useful to post. Do you expect a random difference pattern in fields?

@matzegoebel
Copy link
Contributor Author

ok sounds good. Shall I include the switch and add it to the user guide and tech notes or will someone else do this? What fields are you referring to? I would guess that the differences are random, but I'm not sure.

@dudhia
Copy link
Collaborator

dudhia commented Feb 11, 2021 via email

@matzegoebel
Copy link
Contributor Author

matzegoebel commented Feb 12, 2021 •
edited
Loading

I included a switch called phi_adv_z with default value 1 (original formulation). Do you think that's an appropriate name?
Here are xy plots of phi perturbation and w after 1 hour integration at k=5 (test case em_les again). Looks pretty random but that's hard to say without an in-depth statistical analysis. Let me know if I can do anything else.
w
ph_xy

@dudhia
Copy link
Collaborator

dudhia commented Feb 12, 2021 via email

@matzegoebel
Copy link
Contributor Author

I just figured that I didn't pay attention to get the same initial perturbations for the two em_les runs. Now I manually set the random seed to achieve this. I reproduced the plots that I posted so far. But it doesn't make a big difference. The histogram plots of ph difference get a bit narrower.
hist_ph
scatter_ph
pcolor_ph
hist_all
dmu_timeseries
dmu_timeseries_start
w
ph_xy

@matzegoebel
Copy link
Contributor Author

With such random flows it is difficult, and a more steady flow like the 2d hill case may be more instructive to display this way

Yes, good point. Here are xz difference plots for the em_hill2d_x case (xy plots don't make sense since it's 2D) after 10 hours. The updrafts and downdrafts are about 10% stronger with the new version. This makes sense since the new version employs less averaging for omega. So the modification indeed has some systematic and non-negligible impact.
PH_xz
W_xz

:

@matzegoebel
Copy link
Contributor Author

I also pushed the commit now that includes the switch

@dudhia
Copy link
Collaborator

dudhia commented Feb 15, 2021 via email

@davegill
Copy link
Contributor

@matzegoebel @dudhia
There is a blocky difference at the top of the atmosphere for both ph and w, but more pronounced for ph. Since it is a difference, that points to either the previous or current computation as problematic, correct?

@dudhia
Copy link
Collaborator

dudhia commented Feb 15, 2021 via email

@matzegoebel
Copy link
Contributor Author

The damping options for that case are:
damp_opt = 3
zdamp = 20000. ( ztop = 30000.)
dampcoef = .1

I also thought this looks suspicious. But I guess the damping is too weak to damp away all the waves leading to some strange behavior at the top? But I'm not so familiar with Rayleigh damping to really understand this. If Jimy says it's ok, then I guess it is.

@dudhia
Copy link
Collaborator

dudhia commented Feb 15, 2021 via email

@matzegoebel
Copy link
Contributor Author

matzegoebel commented Feb 16, 2021 •
edited
Loading

yes, I think it has to do with omega=0 at k=kts which introduces a jump if not all waves have been damped away.
I also looked again at the em_les case to check the updraft velocity distributions. The difference between the histogram densities (from minute 30 to 60 for the lowest 20 vertical levels) of the new and old version shows again the stronger up and downdrafts in the new version.
hist_w

@dudhia
Copy link
Collaborator

dudhia commented Feb 18, 2021

This looks fine. The commit message needs to be updated with the addition of the phi_adv_z switch and extra modified files. The new method is 2 and the old one is 1. 1 remains the default.

@matzegoebel
Copy link
Contributor Author

ok I updated the commit message. I also included the last two figures.

@davegill davegill merged commit bd2cde6 into wrf-model:develop Feb 25, 2021
@matzegoebel
Copy link
Contributor Author

@davegill will you add this option to the User Guide?

@dudhia
Copy link
Collaborator

dudhia commented May 5, 2021 via email •
edited
Loading

@weiwangncar
Copy link
Collaborator

@matzegoebel @dudhia We can add this in UG.

vlakshmanan-scala pushed a commit to scala-computing/WRF that referenced this pull request Apr 4, 2024
@matzegoebel @vlakshmanan-scala
Avoid double staggering of ww in vertical advection of geopotential (w…
a1f4f73 
…rf-model#1338)

TYPE: enhancement

KEYWORDS: geopotential,vertical advection

SOURCE: Matthias Göbel (University of Innsbruck)

DESCRIPTION OF CHANGES:
The geopotential equation in WRF is formulated in the following (advective) form where staggering/destaggering is 
marked with overbars:

<img src="https://latex.codecogs.com/svg.latex?\partial_{t}\phi^{\prime}&plus;\mu_{d}^{-1}[m_{x}m_{y}\left(\overline{\overline{U}^\eta\:\delta_{x}\phi}^x&plus;\overline{\overline{V}^\eta\:\delta_{y}\phi}^y\right)&plus;m_{y}\overline{\overline{\Omega}^\eta\:\delta_{\eta}\phi}^\eta-m_{y}gW]=0" title="http://latex.codecogs.com/svg.latex?\partial_{t}\phi^{\prime}+\mu_{d}^{-1}[m_{x}m_{y}\left(\overline{\overline{U}^\eta\:\delta_{x}\phi}^x+\overline{\overline{V}^\eta\:\delta_{y}\phi}^y\right)+m_{y}\overline{\overline{\Omega}^\eta\:\delta_{\eta}\phi}^\eta-m_{y}gW]=0" />

The vertical advection term thus employs a double staggering of Omega. First, Omega is destaggered to mass levels 
and multiplied with the phi gradient. Then the product is staggered back to w-levels. This formulation is analogous to the horizontal terms, but the double averaging of Omega may introduce additional inaccuracies and is not necessary. 

The phi gradient can be first staggered to w-levels and then multiplied with the unaveraged Omega:
 <img src="https://latex.codecogs.com/svg.latex?\overline{\overline{\Omega}^\eta\:\delta_{\eta}\phi}^\eta\rightarrow\Omega\:\overline{\delta_{\eta}\phi}^\eta"  />
For consistency and to avoid numerical instabilities, the same thing has to be done in the acoustic step.

The new namelist option `phi_adv_z` allows switching between the two formulations: 
`phi_adv_z = 1`: old formulation (default)
`phi_adv_z = 2`: new formulation

LIST OF MODIFIED FILES: 
Registry/Registry.EM_COMMON
dyn_em/module_big_step_utilities_em.F
dyn_em/module_small_step_em.F
run/README.namelist 

TESTS CONDUCTED: 
Idealized test cases em_les and em_hill2d_x. 
The vertical velocity is affected by the modification. Updrafts and downdrafts in the em_hill2d_x case become ~10% 
stronger with the new formulation:
![W_xz](https://user-images.githubusercontent.com/17001470/107921009-243d1280-6f6e-11eb-97d0-a897d2d78701.png)
In the LES case stronger updrafts and downdrafts occur more frequently:
![hist_w](https://user-images.githubusercontent.com/17001470/108021629-3c209f00-701f-11eb-80c2-5640201690a0.png)

I came across this issue when trying to transform advective components of temperature, humidity, and momentum 
into Cartesian form in a consistent way, using w for the vertical advection. Without the proposed modifications I 
could not close the budget, because the vertical advection of phi was not consistent with the vertical advection 
in the other budget equations.


RELEASE NOTE: A new namelist option is included that allows a user to avoiding double stagger averaging of Omega in the vertical advection of geopotential. The new namelist option phi_adv_z allows switching between the two formulations: phi_adv_z = 1 (the old formulation, which remains the default), and phi_adv_z = 2 (the new formulation).
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@dudhia dudhia dudhia approved these changes

Assignees
No one assigned
Labels
Develop Branch Dynamics ARW Infrastructure: OK New Feature
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.

None yet
4 participants
@matzegoebel @davegill @dudhia @weiwangncar