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.

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

change edmfx advection #2121

Merged
merged 1 commit into from
Sep 22, 2023
Merged

change edmfx advection #2121

merged 1 commit into from
Sep 22, 2023

Conversation

trontrytel
Copy link
Member

@trontrytel trontrytel commented Sep 20, 2023

Closes #2106

I wasn't sure if we also want to change the tke advection?

@trontrytel trontrytel requested a review from szy21 September 20, 2023 03:54
@trontrytel trontrytel self-assigned this Sep 20, 2023
@trontrytel trontrytel force-pushed the aj/edmfx_advection branch 3 times, most recently from 907a638 to 8cf09db Compare September 20, 2023 04:08
config/model_configs/edmfx_adv_test_box.yml Outdated Show resolved Hide resolved
src/prognostic_equations/advection.jl Outdated Show resolved Hide resolved
toml/edmfx_box_advection.toml Outdated Show resolved Hide resolved
toml/edmfx_box_advection.toml Show resolved Hide resolved
@trontrytel
Copy link
Member Author

trontrytel commented Sep 20, 2023

@tapios , @dennisYatunin - Could you take a look. Below I'm showing results from the edmfx advection test. The first plot shows results for the min_area equal to 1e-2 and the second to 1e-3.

  • The spikes in updraft buoyancy, potential temperature, etc disappear only when min_area is set to 1e-2, which is pretty large.
  • There is a secondary spike in the updraft q_tot (below 5km) that gets larger as min_area gets larger. I don't have an explanation for it right now.
  • The grid mean q_tot doesn't seem to be equal to the area weighted average of the drafts? - Look for example at the point where updraft area is the largest.

TODO - next step: check hat rho^0*q_t^0 + hat rho^1*q_t^1 vs rho q_t on the gm

TODO - next next step: Add the weight for updraft and environment to output. Plot it and see how it looks in problematic areas.

function divide_by_ρa(ρaχ, ρa, ρχ, ρ, turbconv_model)
weight = sgs_weight_function(ρa / ρ, turbconv_model.a_half)
# If ρa = 0, we know that ρa / ρ = 0, which means that weight = 0. However,
# 0 * ρaχ / 0 = NaN, regardless of what ρaχ is, so the linear interpolation
# will always return NaN when ρa = 0. To avoid this problem, we need to add
# a special case for ρa = 0.
return ρa == 0 ? ρχ / ρ : weight * ρaχ / ρa + (1 - weight) * ρχ / ρ
end

min_area = 1e-2

adv_test_2

min_area = 1e-3

adv_test_3

min_area = 1e-4

final_profiles

min_area = 1e-4 Float64
final_profiles

@trontrytel
Copy link
Member Author

Here is another clue. I'm plotting the difference between the grid mean ρq_tot and the area weighted contributions from the drafts computed in two ways: Updraft contribution based on the state Y and based on the cache p.

You can see that the cache based version has a big spike where we have the problematic secondary peak in q_tot.

Updraft computed from the state: (Y.c.ρq_tot - Y.c.sgsʲs[1].ρaq_tot - p.ᶜρa⁰ * p.ᶜspecific⁰.q_tot):

testing_from_Y

Updraft computed from the cache: (Y.c.ρq_tot - Y.c.sgsʲs[1].ρa * p.ᶜspecificʲs[1].q_tot - p.ᶜρa⁰ * p.ᶜspecific⁰.q_tot):

testing_from_p

@trontrytel
Copy link
Member Author

trontrytel commented Sep 21, 2023

Here is how the weight = sgs_weight_function(ρa / ρ, turbconv_model.a_half) looks like for updraft and environment

sgs_weight_function(a, a_half) =
if a < 0
zero(a)
elseif a > 1
one(a)
else
(1 + tanh(2 * atanh(1 - 2 * (1 - a)^(-1 / log2(1 - a_half))))) / 2
end

weight

@tapios
Copy link
Contributor

tapios commented Sep 21, 2023

Here is another clue. I'm plotting the difference between the grid mean ρq_tot and the area weighted contributions from the drafts computed in two ways: Updraft contribution based on the state Y and based on the cache p.

Looks like state and cache are out of sync. May this get better with @charleskawczynski's refactoring of the precomputed quantities to postcomputed quantities?

@szy21
Copy link
Member

szy21 commented Sep 21, 2023

Looks like state and cache are out of sync. May this get better with @charleskawczynski's refactoring of the precomputed quantities to postcomputed quantities?

Maybe. It's more likely from that ρaq_tot is not equal to ρa * q_tot from the regularization I think.

@tapios
Copy link
Contributor

tapios commented Sep 21, 2023

Looks like state and cache are out of sync. May this get better with @charleskawczynski's refactoring of the precomputed quantities to postcomputed quantities?

Maybe. It's more likely from that ρaq_tot is not equal to ρa * q_tot from the regularization I think.

Can we run this test with just the limiters on entrainment/detrainment I sketched, without any additional regularization?

@charleskawczynski
Copy link
Member

charleskawczynski commented Sep 21, 2023

Maybe. It's more likely from that ρaq_tot is not equal to ρa * q_tot from the regularization I think.

Agreed, the cache is currently in sync, the only difference we should see from moving precomputed to post computed quantities is that we should be able to reduce the number of calls by 1 (technically + a fraction due to reducing in callbacks), while keeping the state and cache in sync.

@trontrytel
Copy link
Member Author

Can we run this test with just the limiters on entrainment/detrainment I sketched, without any additional regularization?

This test is only with advection. There is no entrainment/detrainment here

@tapios
Copy link
Contributor

tapios commented Sep 21, 2023

Can we run this test with just the limiters on entrainment/detrainment I sketched, without any additional regularization?

This test is only with advection. There is no entrainment/detrainment here

So entrainment=detrainment=0?

@trontrytel
Copy link
Member Author

So entrainment=detrainment=0?

Yes. The test was meant to look only at advection issues. No other tendencies are applied.

I think that the issue is in the regularization function. I'll try to track it today

@tapios
Copy link
Contributor

tapios commented Sep 21, 2023

So entrainment=detrainment=0?

Yes. The test was meant to look only at advection issues. No other tendencies are applied.

I think that the issue is in the regularization function. I'll try to track it today

Sounds good. Maybe the test is more meaningful with small but nonzero entrainment, detrainment, with the limiter functions on them, to ensure that area fraction remains bounded and updraft properties relax to grid means.

@trontrytel
Copy link
Member Author

Can we merge this change as it is? No matter what we do with the regularization, we still want to switch our advection.

Merging it in would help me with testing for other cases. cc: @tapios @szy21

@szy21
Copy link
Member

szy21 commented Sep 21, 2023

Sounds good to me.

@trontrytel
Copy link
Member Author

bors r+

bors bot added a commit that referenced this pull request Sep 22, 2023
2121: change edmfx advection r=trontrytel a=trontrytel

Closes #2106 

I wasn't sure if we also want to change the tke advection?

2141: Update dependencies r=charleskawczynski a=charleskawczynski

This PR upgrades the dependencies, including ClimaTimeSteppers. The new version fuses some operations, resulting in precision error differences.

Co-authored-by: Anna Jaruga <ajaruga@caltech.edu>
Co-authored-by: Charles Kawczynski <kawczynski.charles@gmail.com>
@bors
Copy link
Contributor

bors bot commented Sep 22, 2023

Build failed (retrying...):

@trontrytel
Copy link
Member Author

bors r+

@bors
Copy link
Contributor

bors bot commented Sep 22, 2023

Already running a review

@bors
Copy link
Contributor

bors bot commented Sep 22, 2023

Build succeeded!

The publicly hosted instance of bors-ng is deprecated and will go away soon.

If you want to self-host your own instance, instructions are here.
For more help, visit the forum.

If you want to switch to GitHub's built-in merge queue, visit their help page.

@bors bors bot merged commit 375443b into main Sep 22, 2023
7 of 8 checks passed
@bors bors bot deleted the aj/edmfx_advection branch September 22, 2023 05:13
@trontrytel
Copy link
Member Author

Below are results with/without using grid mean density Y.c.rho when doing upwind reconstructions for edmfx variables. I'm running the test without any moisture, to highlight that the problem is not with the phase change. The last set of plots is the advection test for our target vertical resolution and using the min_area = 1e3.

  • Changing the density in the advection reconstructions doesn't change the results much.
  • When switching to the target resolution and increasing the min_area, the problem is less severe.

I'm going to merge those settings as the default for prognostic edmfx and move on to the next test case (GABLS).

Using grid mean density in upwind reconstructions:

final_profiles_zoom_in_upwind_recs

Using updraft density in upwind reconstructions:

final_profiles_zoom_in_no_upwind_recs

Target vertical grid:

final_profiles

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Switch advection for edmfx from (rho a) * tracer to rho * (a * tracer)
4 participants