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

Update check_flux_state_consistency() computation #2525

Open
wants to merge 4 commits into
base: master
Choose a base branch
from

Conversation

tcclevenger
Copy link
Contributor

@oksana I update the computation as we talked about.
Notes/Questions:

  1. What is a good way for me to test this? I ran two 5 day simulation with these checks on and off. Neither has mass errors over 1e-10, and both simulations had energy errors around 1e-8, no difference between with or without checks.
  2. Where these checks are performed I don't think we have access to the number of mac_aero_mic steps. The number of steps aren't used anywhere in the shoc code to my knowledge (outside of it being contained in dt). It would be easy to propagate this information, but where are you suggesting we use it?
  3. Should we be using the updated surf_evap_intput in the mass and energy checks for shoc, or the original surf_evap?

@oksanaguba
Copy link
Contributor

tagging here previous discussions #2418 .

@@ -239,7 +240,6 @@ void SHOCMacrophysics::initialize_impl (const RunType run_type)
const auto& pseudo_density = get_field_in("pseudo_density").get_view<const Spack**>();
const auto& omega = get_field_in("omega").get_view<const Spack**>();
const auto& surf_sens_flux = get_field_in("surf_sens_flux").get_view<const Real*>();
const auto& surf_evap = get_field_in("surf_evap").get_view<const Real*>();
Copy link
Contributor

Choose a reason for hiding this comment

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

needs to be grabbed back

@@ -376,7 +377,7 @@ void SHOCMacrophysics::initialize_impl (const RunType run_type)
const auto& water_flux = get_field_out("water_flux").get_view<Real*>();
const auto& ice_flux = get_field_out("ice_flux").get_view<Real*>();
const auto& heat_flux = get_field_out("heat_flux").get_view<Real*>();
shoc_postprocess.set_mass_and_energy_fluxes (surf_evap, surf_sens_flux,
shoc_postprocess.set_mass_and_energy_fluxes (surf_evap_input, surf_sens_flux,
Copy link
Contributor

Choose a reason for hiding this comment

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

use original surf_evap here.

@@ -499,17 +506,16 @@ void SHOCMacrophysics::check_flux_state_consistency(const double dt)
{
using PC = scream::physics::Constants<Real>;
const Real gravit = PC::gravit;
const Real qmin = 1e-12; // minimum permitted constituent concentration (kg/kg)
Copy link
Contributor

Choose a reason for hiding this comment

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

@bogensch @AaronDonahue @PeterCaldwell this lower bound was taken from EAM setup (iirc EAM won't run if this bount is set to 0.0). Do we need a certain nontrivial bound in scream or can we always assume lower acceptable bounds for tracers are zero. In this case it is only for vapor. If a certain bound is needed, Conrad will put it in physics constants.

Copy link
Contributor

Choose a reason for hiding this comment

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

maybe it is a question for @brhillman instead, since in EAM the comment mentions radiation

    real(r8),intent(in)    :: qminc  ! minimum value of mass mixing ratio (kg/kg)
                                     ! normally 0., except water 1.E-12, for radiation.

const auto adjust = cc*qv_i(k)*pseudo_density_i(k)/mm;
qv_i(k) = (qv_i(k)*pseudo_density_i(k) - adjust)/pseudo_density_i(k);
const auto excess_mass = surf_evap_mass + column_vapor_mass(last_pack_idx)[last_pack_entry];
if (excess_mass < 0) {
Copy link
Contributor

Choose a reason for hiding this comment

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

qmin?

// reciprocal of pseudo_density at the bottom layer
const auto rpdel = 1.0/pseudo_density_i(last_pack_idx)[last_pack_entry];
// Define lambda to calculate column vapor mass
auto column_vapor_mass = [&](const int k) {
Copy link
Contributor

Choose a reason for hiding this comment

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

add function in the mname then?

"than mass of surf_evap.\n");

// Redistribute excess mass
Kokkos::parallel_for(Kokkos::TeamVectorRange(team, nlev_m_1_packs), [&](const int& k) {
Copy link
Contributor

Choose a reason for hiding this comment

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

excess mass needs to be used here

@oksanaguba
Copy link
Contributor

@tcclevenger here is the code i ran in EAM (in place of qneg4 call in physpkg line 1767)

!vapor hard set bound
bb = 1e-12_r8
!convert to/from Pa
ffactor = ztodt*gravit

do i=1,ncol
  !convert cflx in Pa units
  cc = cam_in%cflx(i,1)*ffactor

  qqdp1 = state%q(i,pver,1)*state%pdel(i,pver)
  !if adding cflx to the bottom violates the bound
  if ( (cc + qqdp1) < bb) then

    !print *, 'OG cc value is ', cc
    cc=abs(cc)

    !comute total Pa mass of water vapor in column
    mm=0.0_r8
    do k=1,pver
      mm = mm + state%q(i,k,1)*state%pdel(i,k)
    enddo
    mass_before = mm + cam_in%cflx(i,1)*ffactor

    !set cflx to the value that won't violate bb bound
    cam_in%cflx(i,1) = (bb - qqdp1)/ffactor

    if (mm < cc) then
      !abort, cflx value is negative and by abs is larger than total mass of vapor in column
      stop
    else
      !compute mass ot 1:plev-1 column
      mm = mm - qqdp1
      !redistribute excess = -cc - (bb - qqdp1)
      !check: excess + cflx*ffactor = -cc (this is what we want)
      !cc is negative excess now
      cc = -cc - (bb - qqdp1)
      do k=1,pver-1
        qq = state%q(i,k,1)
        pp = state%pdel(i,k)
        adjust = cc *  qq*pp/mm
        state%q(i,k,1) = (qq*pp + adjust)/pp
      enddo
    endif
!print *, "OG cc, mm, new mm",cc,mm,sum(state%q(i,:,1)*state%pdel(i,:))
!print *, 'OG mass before, after',mass_before, (sum(state%q(i,:,1)*state%pdel(i,:)) + (bb - qqdp1))
!print *, 'OG mass before - after',mass_before - (sum(state%q(i,:,1)*state%pdel(i,:)) + (bb - qqdp1))
!with right code rel error is ~1e-16
!print *, 'OG rel mass before - after',(mass_before - (sum(state%q(i,:,1)*state%pdel(i,:)) + (bb - qqdp1)))/mass_before
  endif
enddo

@oksanaguba
Copy link
Contributor

oksanaguba commented Sep 29, 2023

The code sets cflx to an acceptable value and redistributes leftovers in the column 1:plev-1. The mass check for the code is rel error ~1e-16.
Diagnostics for qneg4 and the new code is in this plot

Screen Shot 2023-09-29 at 3 27 36 PM

Here green line is the old code, blue line is the new. Diagnostics is dWater minus (expected fluxes). Fluxes values for these runs are ~1e-3 to 1e-1. One interesting detail is that the run above is for macmic=1, not for the default macmic=6. With macmic=6, the diagnostics for the new code is actually worse. In EAM case the code does not take into account that there are many iterations of macmic loop; i did not try to understand it better than that, i just switched to 1.

@oksanaguba
Copy link
Contributor

i forgot about mass_borrower for the runs above, will repeat them.

@AaronDonahue
Copy link
Contributor

@oksanaguba and @tcclevenger , what is the status of this PR? Is it still WIP?

@tcclevenger
Copy link
Contributor Author

@AaronDonahue For now still WIP, I need to implement Oksana's suggestions. I'll try and get to that this week.

@E3SM-Bot
Copy link
Collaborator

Status Flag 'Pre-Test Inspection' - Auto Inspected - Inspection is Not Necessary for this Pull Request.

@E3SM-Bot
Copy link
Collaborator

Status Flag 'Pull Request AutoTester' - Failure: Timed out waiting for job SCREAM_PullRequest_Autotester_Mappy to start: Total Wait = 1803

  • Other jobs have been previously started - We must stop them...

@E3SM-Bot
Copy link
Collaborator

Status Flag 'Pull Request AutoTester' - Failure: Timed out waiting for job SCREAM_PullRequest_Autotester_Weaver to start: Total Wait = 1803

  • Other jobs have been previously started - We must stop them...

@tcclevenger tcclevenger linked an issue Oct 12, 2023 that may be closed by this pull request
@tcclevenger
Copy link
Contributor Author

@oksanaguba Could you take a look at this again and maybe test? I updated the computation based on your code.

@oksanaguba
Copy link
Contributor

@tcclevenger thanks, i cannot do it immediately, maybe in 1-2 weeks.

@AaronDonahue
Copy link
Contributor

@tcclevenger , is there any reason to keep this PR open? If so, can it be revived and updated w/ the latest master?

@tcclevenger
Copy link
Contributor Author

@tcclevenger , is there any reason to keep this PR open? If so, can it be revived and updated w/ the latest master?

@oksanaguba Is there any reason to keep this open? Can't remember if we were going a different direction.

@oksanaguba
Copy link
Contributor

@AaronDonahue i haven't looked at this for a while. i remember that it required nontrivial changes, because we want to use limiter on cflx each macmic iteration, not just once before macmic loop.

i will try to talk to Conrad soon to figure this out.

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.

ENERGY, keep open: Convert qneg4 routine to C++
4 participants