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

Specified Monin-Obukhov Length #1418

Open
wants to merge 35 commits into
base: main
Choose a base branch
from
Open

Conversation

hgopalan
Copy link
Contributor

@hgopalan hgopalan commented Dec 20, 2024

There are three ways to specify non-neutral stratification for ABL simulations: (i) specified heat-flux, (ii) specified temperature or heating/cooling rate and (iii) specified Monin-Obukhov length. This PR adds the capability for the third methodology in which the Monin-Obukhov length is kept constant. This method is used often in RANS simulations to do wind turbine siting for non-neutral cases (faster convergence rate) and minimize the mismatch between wall boundary condition and inflow profile.

This method is also robust for running SBL cases in which the non-linear iteration for MOL may not converge.

Please check the type of change introduced:

  • Bugfix
  • [X ] Feature
  • Code style update (formatting, renaming)
  • Refactoring (no functional changes, no api changes)
  • Build related changes
  • Documentation content changes
  • Other (please describe):

Checklist

The following is included:

  • new unit-test(s)
  • new regression test(s)
  • documentation for new capability

This PR was tested by running:

  • the unit tests
    • on GPU
    • [X ] on CPU
  • the regression tests
    • on GPU
    • on CPU

@hgopalan hgopalan added the enhancement New feature or request label Dec 20, 2024
@hgopalan hgopalan marked this pull request as ready for review February 24, 2025 05:38
@hgopalan
Copy link
Contributor Author

Need additional test cases for immersed forcing

amrex::Real m_gamma_h{5.0};
amrex::Real m_beta_h{16.0};

static amrex::Real stability(
Copy link
Contributor

Choose a reason for hiding this comment

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

This function seems to be getting repeated implementations across multiple files? Just have one and call those?

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 did not know how to do that across namespaces so created one at every place.

Copy link
Contributor

Choose a reason for hiding this comment

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

You could just put it in a common place and call it across the places you need. For instance you could add these in some type of wall function utilities. We don't want to be repeating code.

Copy link
Contributor

Choose a reason for hiding this comment

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

Is the rans_1d.info file identical across the cases? You might be able to just have one and reuse it for the other tests?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

No they are different for each case.

const amrex::Real kappa = 0.41;
const amrex::Real cd_max = 1000.0;
amrex::Real non_neutral_neighbour = 0.0;
Copy link
Contributor

Choose a reason for hiding this comment

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

Can you const these with a ternary expression?

uTarget * ux2r / (tiny + std::sqrt(ux2r * ux2r + uy2r * uy2r));
const amrex::Real uyTarget =
uTarget * uy2r / (tiny + std::sqrt(ux2r * ux2r + uy2r * uy2r));
ustar / kappa * (std::log(0.5 * dx[2] / z0) - non_neutral_cell);
Copy link
Contributor

Choose a reason for hiding this comment

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

Are these calcs/the model going to be in the documentation? Not sure how to interpret non_neutral_cell and non_neutral_neighbor.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The walk through talks about the details already. The non_neutral are stability function corrections.

const auto& geom = m_mesh.Geom(lev);
const auto& dx = geom.CellSizeArray();
const amrex::Real drag_coefficient = m_drag_coefficient / dx[2];
amrex::FArrayBox ref_theta_fab(bx, 1, amrex::The_Async_Arena());
amrex::Array4<amrex::Real> const& ref_theta_arr = ref_theta_fab.array();
m_transport.ref_theta_impl(lev, mfi, bx, ref_theta_arr);
const amrex::Real gravity_mod = 9.81;
Copy link
Contributor

Choose a reason for hiding this comment

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

Do you want to be using the value of gravity from the input file? We have incflo.gravity or some such right?

//! We do not know the actual temperature so use cell above
const amrex::Real thetastar =
theta * ustar * ustar / (kappa * gravity_mod * mol_length);
const amrex::Real surf_temp =
Copy link
Contributor

Choose a reason for hiding this comment

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

Is this some type of interpolation to get to surf temp? Wondering why this is using theta2 vs say theta or the temperature at the bottom cell?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

For terrain, the value at the point above the cell is computed from the solution as it is not a wall. However, we need to force a MOST boundary condition as a forcing term. So we correct the theta at the cell using the cell above it by computing what is the MOST solution. We already had this in place for wind speed. It gets more complicated with temperature.

amrex::Real m_gamma_m{5.0};
amrex::Real m_beta_m{16.0};

static amrex::Real stability(
Copy link
Contributor

Choose a reason for hiding this comment

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

Same comment as repeated definition of this function

if (zeta > 0) {
psi_m = -m_gamma_m * zeta;
} else {
amrex::Real x =
Copy link
Contributor

Choose a reason for hiding this comment

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

This variable name is confusing, this is not plo + i*dx.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Copy link
Contributor

Choose a reason for hiding this comment

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

Why not just call that function?

Copy link
Contributor

Choose a reason for hiding this comment

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

Also context matters, that is a standalone small function where x isn't going to get confused with the spatial direction. Which isn't the case in your code.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Why not just call that function?

I want to have the flexibility of using a different stability function in the future within the context of this model. Currently it is set to the same as the one in MODData.H.

Copy link
Contributor

Choose a reason for hiding this comment

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

In that case, create a function for calculating psi_m/psi_h in your namespace and call the MOData function within that. Or simply just call the MOData function for now, and replace it when you add a different stability function.

const amrex::Real vv = vold_arr(i, j, k, 1);
const amrex::Real wspd =
std::sqrt(uu * uu + vv * vv);
const amrex::Real theta2 = told_arr(i, j, k);
Copy link
Contributor

Choose a reason for hiding this comment

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

maybe a more descriptive variable name? you have theta2 to refer to theta at k+1 in a different place.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

In this case the value at k is correct so I am referring to it as theta2.

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

Successfully merging this pull request may close these issues.

2 participants