Skip to content
This repository has been archived by the owner on Mar 1, 2023. It is now read-only.

Diagnostics module #508

Merged
merged 1 commit into from
Dec 17, 2019
Merged

Diagnostics module #508

merged 1 commit into from
Dec 17, 2019

Conversation

kpamnany
Copy link
Contributor

Refactored implementation of the diagnostics prototyped in #467.

This introduces the first version of an infrastructure for gathering diagnostics on the DG grid and moves the diagnostics implemented in #467 onto that infrastructure. It also moves the code into a CLIMA module.

Diagnostics gathering is still triggered by the callback mechanism (and must be registered in the driver) until the proper timestepping mechanism is implemented.

Output is to JLD2 and VariableTemplates can be used with the help of the *_vars.jl files introduced here for ease in post-processing.

In future PRs:

  • More abstraction and integration with Oceananigans diagnostics
  • Reducing MPI communication

@slifer50 and @smarras79: please check the output for correctness relative to the previous PR.

Cc: @tapios, @sandreza, @christophernhill, @blallen, @glwagner

@kpamnany
Copy link
Contributor Author

@jkozdon's comment that the "sums should be integral approximations (involving quadrature weights and jacobian determinants)" still applies. That is hopefully easier to implement now.

@codecov
Copy link

codecov bot commented Oct 30, 2019

Codecov Report

Merging #508 into master will increase coverage by 3.13%.
The diff coverage is 83.17%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #508      +/-   ##
==========================================
+ Coverage   68.02%   71.15%   +3.13%     
==========================================
  Files          82       83       +1     
  Lines        5620     5814     +194     
==========================================
+ Hits         3823     4137     +314     
+ Misses       1797     1677     -120
Impacted Files Coverage Δ
src/CLIMA.jl 100% <ø> (ø) ⬆️
src/Diagnostics/graph_diagnostic.jl 0% <0%> (ø)
src/Diagnostics/diagnostic_vars.jl 100% <100%> (ø)
...c/Utilities/VariableTemplates/VariableTemplates.jl 94.87% <80%> (+0.08%) ⬆️
src/Diagnostics/Diagnostics.jl 91.52% <91.52%> (ø)
src/DGmethods/DGmodel_kernels.jl 97.14% <0%> (+0.19%) ⬆️
src/Arrays/MPIStateArrays.jl 68.87% <0%> (+6.12%) ⬆️
... and 8 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a1bc4b2...d0c81ec. Read the comment docs.

@smarras79
Copy link
Contributor

@slifer50 and @smarras79: please check the output for correctness relative to the previous PR.

will try it between later today and tomorrow and let you. Thanks for making this a module!

H = Rm_sfc * T_sfc / grav;
p = P_sfc * exp(-xvert/H);
#Density, Temperature
TS = LiquidIcePotTempSHumEquil_no_ρ(θ_liq, q_pt, p)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will need to be updated after PR #506 is merged. @charleskawczynski can help.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it possible to get density here and pass it in? Or do we need a LiquidIcePotTempSHumEquil_given_pressure(θ_liq, q_pt, p) constructor? I tried looking into this recently and ran into some problems with the equations.

Copy link
Contributor

@smarras79 smarras79 Nov 5, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you can get T from theta_l and q_t using your functions then you can get rho from the equation of state using p and T as I do below in my matlab script (@tapios please, confirm that the operations below would be correct here):

  `p = p_sfc * exp(-z / H)

%Exner
exner = (p / p_sfc)^(R_d/cp_d);

%T, Tv 
T     = exner*theta_liq + Lv0 * q_liq / (cpm*exner);

%Density
rho  = p/(Rm*T);
`

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you can get T from theta_l and q_t using your functions then you can get rho from the equation of state using p and T as I do below in my matlab script (@tapios please, confirm that the operations below would be correct here):

  `p = p_sfc * exp(-z / H)

%Exner
exner = (p / p_sfc)^(R_d/cp_d);

This looks ok.

%T, Tv
T = exnertheta_liq + Lv0 * q_liq / (cpmexner);

The second term should not be divided by the Exner function: Just

T = exner*theta_liq + Lv0 * q_liq / (cpm);

(I put the derivations in the design docs; please refer to those.)

@charleskawczynski I think providing a function LiquidIcePotTempSHumEquil_given_pressure(θ_liq, q_pt, p) would be useful. In essence, this is how @smarras79 is calculating T.

%Density
rho = p/(Rm*T);
`

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just for reference, we'll have to use LiquidIcePotTempSHumNonEquil_given_pressure for now until we come up with a consistent constructor for an equilibrium phase.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This driver is now outdated and needs to be replaced with the currently running one.

@smarras79
Copy link
Contributor

smarras79 commented Nov 5, 2019 via email

@kpamnany kpamnany mentioned this pull request Nov 5, 2019

e_int = e̅_tot - 1//2 * (u̅^2 + v̅^2 + w̅^2) - grav * zvals[k,ev]

ts = PhaseEquil(convert(FT, e_int), q̅_tot, state.ρ)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might be nice if we could pass in moist::MoistureModel, orientation::Orientation, state::Vars, aux::Vars and dispatch to the dry or moist thermo state functions instead, but maybe that'll be the next step.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agree! We have state and can easily have aux. Would moist be this and orientation be this? Or would they need to come from somewhere else?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yep, those are the correct structs.

Comment on lines 185 to 191
u̅ = state.ρu[1] / state.ρ
v̅ = state.ρu[2] / state.ρ
w̅ = state.ρu[3] / state.ρ
q̅_tot = state.moisture.ρq_tot / state.ρ
ũ = ha.ρu / ha.ρ
ṽ = ha.ρv / ha.ρ
w̃ = ha.ρw / ha.ρ
ẽ = ha.e_tot / ha.ρ
q̃_tot = ha.ρq_tot / ha.ρ

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would suggest changing the notation to avoid confusion. The tilde and overbar look too similar and can lead to coding mistakes.
For example:
mean u --> \underline{u}
instantaneous u: --> u
fluctuation of u: --> \tilde{u}

or anything that clearly distinguishes each quantity

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't have a preference -- your thoughts @akshaysridhar?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agree with instantaneous u: --> u , this affects 185 to 188 and quantities dependent on these variables. Fluctuating quantities seem fine with \prime to me. No preference for the rest.

@kpamnany kpamnany force-pushed the kp/diagnostics branch 2 times, most recently from baebf28 to 1e2d9cd Compare November 18, 2019 22:45
Comment on lines 14 to 41
z::FT
# state and functions of state
u::FT
v::FT
w::FT
q_tot::FT
e_tot::FT
q_liq::FT
θ::FT
θ_liq::FT
θ_v::FT
e_int::FT
# vertical fluxes
w′ρ′::FT
w′u′::FT
w′v′::FT
w′q_tot′::FT
w′q_liq′::FT
w′q_vap′::FT
w′θ′::FT
w′θ_v′::FT
w′θ_liq′::FT
# variances
u′u′::FT
v′v′::FT
w′w′::FT
# skewness
w′w′w′::FT
Copy link
Contributor

@smarras79 smarras79 Nov 19, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Although we use symbols within the main clima code, it would be better to not use symbols for the variables that are written to file.
Not every user may be able to use symbols when calling @kpamnany 's new graph_diagnostics.jl from the command line
(Notice the graph_diagnostics.jl is passed the list of variables to be plotted:
>> julia graph_diagnostics.jl <diagnostic_file.jld2> <diagnostic_name> where
is a list of the variables to be plotted.

In atmospheric software packages we usually find a combination of short and sometime full names such as:
"TH", "Potential temperature"
"THV", "Virtual potential temperature"
"PRS", "Pressure"
"rho", "Air density"
"qv", "Water vapor mixing ration"
etc. (these are examples from CM1 by Bryan, but they all more or less use a similar convention)
Bottom line, I would totally avoid greek symbols to name the diagnostic variables.

Suggestions (@kpamnany @akshaysridhar @tapios):

<w′θ′> --> vert_eddy_heat_flx
<wθ> --> vert_heat_flx
<w′q_tot>′ --> vert_eddy_qt_flx
<wq_tot> --> vert_qt_flx
<u′u′> --> uvariance
<u′u′> --> vvariance
<w′w′> --> wvariance
<w′u′> --> vert_eddy_u_flx
<w′u′> --> vert_eddy_v_flx

etc.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added altvars_diagnostic, as a mirror of vars_diagnostic, with some of the changes suggested here. If you give me a full list of what you'd like changed, I'll make the changes and then update graph_diagnostic.jl to use altvars_diagnostic instead, which will allow you to use non-unicode names in post-processing while still allowing the CLIMA code to use unicode.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not to make this go back and forth but, I think a better solution here would be to write mapping functions for each symbol/expression, and calls to print/show can be intercepted by the mapping funcs, instead of having two structs of the same thing.

Copy link
Member

@charleskawczynski charleskawczynski Nov 22, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or use something like:

struct VerboseName{T} end
VerboseName(T) = VerboseName{T}()

nice_name(::VerboseName{:θ})        = :theta
nice_name(::VerboseName{:e_int})    = :e_int

struct Foo{FT}
  θ::FT
  e_int::FT
end

nice_names(f::Foo) = nice_name.(VerboseName.(fieldnames(typeof(f))))

f = Foo(1.0,2.0)

@show nice_names(f)

Copy link
Contributor

@tapios tapios Nov 22, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed. In the diagnostics issue, I outlined what we’ll need in the end: a table that defines short and verbose names for the output fields, and units, to be able to generate self-documenting output. Short names will be variable names in files. May not need to be now, but that’s where we’ll need to go eventually.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Having primes etc. in the code is fine for now though. I wouldn’t go through the trouble to create more conventional variable names in the code right now.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll just switch to the conventional variable names now, since we want to do that anyway.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll just switch to the conventional variable names now, since we want to do that anyway.

If by conventional you mean the long names @smarras79 suggested, my take is that is not needed right now. There are no general conventions for the names, and your primed notation had the advantage of clarity for now. Most importantly, it think it would be good to get this merged, and then we can refine.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Simone has been working with a modified version of this code with plain text names due to difficulties with specifying the Unicode names we're using. I've just switched to his version; not a difficult thing. I agree on getting this merged ASAP!

@kpamnany
Copy link
Contributor Author

@slifer50 has written a test for this PR per the suggestion from @tapios but it doesn't work too well in the absence of #545. So we're waiting on that to be merged after which we will update this and also add the test.

@kpamnany kpamnany force-pushed the kp/diagnostics branch 6 times, most recently from 01b3fd7 to 26053a1 Compare December 13, 2019 21:28
@kpamnany
Copy link
Contributor Author

@simonbyrne: this is good to go.

@simonbyrne
Copy link
Member

bors r+

bors bot added a commit that referenced this pull request Dec 16, 2019
508: Diagnostics module r=simonbyrne a=kpamnany

Refactored implementation of the diagnostics prototyped in #467.

This introduces the first version of an infrastructure for gathering diagnostics on the DG grid and moves the diagnostics implemented in #467 onto that infrastructure. It also moves the code into a CLIMA module.

Diagnostics gathering is still triggered by the callback mechanism (and must be registered in the driver) until the proper timestepping mechanism is implemented.

Output is to JLD2 and `VariableTemplates` can be used with the help of the `*_vars.jl` files introduced here for ease in post-processing.

In future PRs:
- More abstraction and integration with Oceananigans diagnostics
- Reducing MPI communication

@slifer50 and @smarras79: please check the output for correctness relative to the previous PR.

Cc: @tapios, @sandreza, @christophernhill, @blallen, @glwagner 


Co-authored-by: Kiran Pamnany <kpamnany@gmail.com>
@bors
Copy link
Contributor

bors bot commented Dec 16, 2019

Build failed

@kpamnany
Copy link
Contributor Author

bors r+

bors bot added a commit that referenced this pull request Dec 16, 2019
508: Diagnostics module r=kpamnany a=kpamnany

Refactored implementation of the diagnostics prototyped in #467.

This introduces the first version of an infrastructure for gathering diagnostics on the DG grid and moves the diagnostics implemented in #467 onto that infrastructure. It also moves the code into a CLIMA module.

Diagnostics gathering is still triggered by the callback mechanism (and must be registered in the driver) until the proper timestepping mechanism is implemented.

Output is to JLD2 and `VariableTemplates` can be used with the help of the `*_vars.jl` files introduced here for ease in post-processing.

In future PRs:
- More abstraction and integration with Oceananigans diagnostics
- Reducing MPI communication

@slifer50 and @smarras79: please check the output for correctness relative to the previous PR.

Cc: @tapios, @sandreza, @christophernhill, @blallen, @glwagner 


Co-authored-by: Kiran Pamnany <kpamnany@gmail.com>
@bors
Copy link
Contributor

bors bot commented Dec 16, 2019

Build failed

… the state array as well as some notable second moments average along the horizontal plane for each vertical level

Refactored into CLIMA.Diagnostics. Add template for post-processing and
graphing diagnostics output. Add test case. Support multiple ranks.
@kpamnany
Copy link
Contributor Author

bors r+

bors bot added a commit that referenced this pull request Dec 17, 2019
508: Diagnostics module r=kpamnany a=kpamnany

Refactored implementation of the diagnostics prototyped in #467.

This introduces the first version of an infrastructure for gathering diagnostics on the DG grid and moves the diagnostics implemented in #467 onto that infrastructure. It also moves the code into a CLIMA module.

Diagnostics gathering is still triggered by the callback mechanism (and must be registered in the driver) until the proper timestepping mechanism is implemented.

Output is to JLD2 and `VariableTemplates` can be used with the help of the `*_vars.jl` files introduced here for ease in post-processing.

In future PRs:
- More abstraction and integration with Oceananigans diagnostics
- Reducing MPI communication

@slifer50 and @smarras79: please check the output for correctness relative to the previous PR.

Cc: @tapios, @sandreza, @christophernhill, @blallen, @glwagner 


Co-authored-by: Kiran Pamnany <kpamnany@gmail.com>
@bors
Copy link
Contributor

bors bot commented Dec 17, 2019

Build succeeded

@bors bors bot merged commit d0c81ec into master Dec 17, 2019
@bors bors bot deleted the kp/diagnostics branch December 17, 2019 23:54
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants