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

Machine Learning based vertical diffusivity in EPBL mixing scheme used for ocean surface boundary layer #737

Open
wants to merge 40 commits into
base: dev/gfdl
Choose a base branch
from

Conversation

aakashsane
Copy link

@aakashsane aakashsane commented Oct 9, 2024

Requesting to add a few subroutines to the EPBL vertical mixing module. 
These changes enhances the existing vertical diffusivity used in EPBL with machine learning, and makes it constrained on second moment closure. Using symbolic regression and empirical fitting, a shape function ( g(\sigma) ) has been formulated that responds to changes in the surface forcing (Latitude, wind stress, surface buoyancy flux, boundary layer depth). g(\sigma) goes from zero to 1 and its skewness changes as per surface forcing conditions. The velocity scale, v_0, is an approximation that depends on (Coriolis f, ustar, and surface buoyancy flux).
When v_0 is combined with g(\sigma) and multiplied by the energetics based boundary layer depth h, i.e \nu = . g(\sigma) X v_0 X h, we get a diffusivity which is constrained on a second moment closure.

The subroutines are activated by using the flags:

  1. Equation_Discovery_shape = True
  2. Equation_Discovery_velocity = True
  3. Equation_Discovery_velocity_h = False
  4. activates the new equation for shape function. 2. Activates velocity scale equations that approximate neural networks as described in Sane et al. 2023. 3. Is a new velocity scale that includes boundary layer depth as input and agrees with well known scaling.
    Either 2. or 3. should be 'True', both cannot be True or False.

The new subroutines have been tested by running the OM4 configuration (Adcroft et al. 2019) and comparing the simulation against observations (SST from WOA and Mixed Layer Depth from ARGO). The equations are approximating the neural networks. The neural network based diffusivity has been published in Sane et al. 2023 ( https://doi.org/10.1029/2023MS003890 ) and the equations based diffusivity will soon be submitted for a publication.
The changes have also been tested using the scale test.

All the commits can be squashed into one as only the latest is relevant.

aakashsane and others added 30 commits July 17, 2024 22:37
updating ML_diffusivity to latest dev/gfdl
replaced the sigma to z coord algorithm used to map neural network output with a simpler and better algorithm that finds the shape function on the hz vertical grid.
some typo corrections
f_lower is a lower cap on abs_f used inside equation for v_0. 
A cap is required to avoid singularity. Capping at any value below 1 deg is okay, solution is not sensitive. 
f_lower was tested for 1 deg and 0.1 deg. 
SST and MLD did not change.
updating latest dev/gfdl with ML_diffusivity
Copy link

codecov bot commented Oct 31, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 40.94%. Comparing base (f90b071) to head (992f266).
Report is 5 commits behind head on dev/gfdl.

Additional details and impacted files
@@             Coverage Diff              @@
##           dev/gfdl     #737      +/-   ##
============================================
+ Coverage     36.63%   40.94%   +4.30%     
============================================
  Files           274       42     -232     
  Lines         84153     5288   -78865     
  Branches      15834     1013   -14821     
============================================
- Hits          30829     2165   -28664     
+ Misses        47509     2938   -44571     
+ Partials       5815      185    -5630     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

Copy link
Member

@adcroft adcroft left a comment

Choose a reason for hiding this comment

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

I can't (shouldn't) make the final approval on this since I'm involved in the research, but there are some obvious things that need to be changed before someone else reviews the PR. Once these are cleaned up, my take is it will be smooth sailing.

You added several files that I think should not be included:

  • .testing/tc4/Makefile is generated by autoconf using the Makefile.in template
  • .testing/gen_data.dSYM/Contents.Info.plist - looks like something leftover from your a mac build process
  • .testing/gen_grid.dSYM/Contents.Info.plist - ditto

Please remove these files. I believe the only file you meant to submit changes for is MOM_energetic_PBL.F90. There's also a real variable with missing unit documentation which I added a separate inline comment for.



! variables for ML based diffusivity
real :: v0_dummy ! Variable which get recycled, set equal to v_0
Copy link
Member

Choose a reason for hiding this comment

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

This real variable needs units in the comment (looks like it should be the same as CS%vo)

Copy link
Author

@aakashsane aakashsane Nov 21, 2024

Choose a reason for hiding this comment

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

Thanks for pointing it out.
I have completed:

  1. Units added to that variable.
  2. All the files mentioned above have been deleted.

@adcroft
Copy link
Member

adcroft commented Nov 25, 2024

I've just been discussing with Bob, and we think we need to add a new version of get_param to solve a problem that has been highlighted by this PR, namely that we're adding 16 run-time parameters, for one scheme, when they could be a single line. We'll look into adding a vector default option and get back to this PR shortly.

if (CS%MixLenExponent==2.0) then
MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * &
(max(0.0, (MLD_guess - dz_rsum)*I_MLD) )**2 ! CS%MixLenExponent

Copy link
Member

Choose a reason for hiding this comment

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

It looks like you are replacing some existing code, rather than adding a new option. I suspect you need to add elseif (CS%eqdisc_v0 .or.CS%eqdisc_v0h) then

Copy link
Author

Choose a reason for hiding this comment

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

The subroutines associated with CS%eqdisc_v0 and CS%eqdisc_v0h are different.
CS%eqdisc_v0 calls the subroutine 'get_eqdisc_v0(CS,absf,B_flux,u_star,v0_dummy)'
while
CS%eqdisc_v0h calls 'get_eqdisc_v0h(CS,B_flux,u_star,MLD_guess,v0_dummy)'.

Both of them formulate v0 in a different manner.
So it would not be ideal to call them in a single elseif statement.

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

Successfully merging this pull request may close these issues.

3 participants