Skip to content

Commit

Permalink
Nudging Mods (#1492)
Browse files Browse the repository at this point in the history 
* Best nudging FLOWMAS simulation yet.

* Clean up and add doc description.

* Clarify docs.

---------

Co-authored-by: Aaron Lattanzi <[email protected]>
Co-authored-by: Ann Almgren <[email protected]>
  • Loading branch information
@AMLattanzi @asalmgren
3 people authored Mar 13, 2024
1 parent e56589c commit 2e777a2
Show file tree
Hide file tree
Showing 5 changed files with 39 additions and 32 deletions.
34 changes: 18 additions & 16 deletions Docs/sphinx_doc/BoundaryConditions.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,17 +139,22 @@ It is important to note that external Dirichlet boundary data should be specifie
as the value on the face of the cell bounding the domain, even for cell-centered
state data.

Specified Domain BCs
Real Domain BCs
----------------------

When we use specified lateral boundary conditions, we read time-dependent Dirichlet data
from a file. The user may specify (in the inputs file)
When using real lateral boundary conditions, time-dependent observation data is read
from a file. The observation data is utilized to directly set Dirichlet values on the
lateral domain BCs as well as nudge the solution state towards the observation data.
The user may specify (in the inputs file)
the total width of the interior Dirichlet and relaxation region with
``erf.real_width = <Int>`` (yellow + blue)
and analogously the width of the interior Dirichlet region may be specified with
``erf.real_set_width = <Int>`` (yellow).
We note that all dycore variables are set and relaxed, but moisture and other scalars
are only set in the yellow region if present in the boundary file.
The real BCs are only imposed for :math:`\psi = \left{ \rho; \; \rho \theta; \; \rho q_v; \; u; \; v \right}`.
Due to the staggering of scalars (cell center) and velocities (face center) with an Arakawa C grid,
we reduce the relaxation width of the scalars :math:`\left{ \rho; \; \rho \theta; \; \rho q_v \right}` by 1
to ensure the momentum updates at the last relaxation cell involve a pressure gradient that is computed with
relaxed and non-relaxed data.

.. |wrfbdy| image:: figures/wrfbdy_BCs.png
:width: 600
Expand All @@ -166,26 +171,23 @@ are only set in the yellow region if present in the boundary file.

.. _`Skamarock et al. (2021)`: http://dx.doi.org/10.5065/1dfh-6p97

Here we describe the relaxation algorithm.

Within the interior Dirichlet region (yellow), the RHS is exactly 0.
However, within the relaxation region (blue), the RHS (:math:`F`) is given by the following:
Within the interior Dirichlet cells, the RHS is exactly :math:`\psi^{n} - \ps^{BDY} / \Delta t`
and, as such, we directly impose this value in the yellow region.
Within the relaxation region (blue), the RHS (:math:`F`) is given by the following:

.. math::
\begin{align}
F &= G + R, \\
\psi^{\prime} &= \psi^{n} + \Delta t \; G, \\
R &= H_{1} \left( \psi^{FP} - \psi^{\prime} \right) - H_{2} \Delta^2 \left( \psi^{FP} - \psi^{\prime} \right), \\
R &= \left[ H_{1} \left( \psi^{BDY} - \psi^{\*} \right) - H_{2} \Delta^2 \left( \psi^{BDY} - \psi^{\*} \right) \right] \exp \left(-C_{01} \left(n - {\rm SpecWidth}\right) \right), \\
H_{1} &= \frac{1}{10 \Delta t} \frac{{\rm SpecWidth} + {\rm RelaxWidth} - n}{{\rm RelaxWidth} - 1}, \\
H_{2} &= \frac{1}{50 \Delta t} \frac{{\rm SpecWidth} + {\rm RelaxWidth} - n}{{\rm RelaxWidth} - 1},
\end{align}
where :math:`G` is the RHS of the evolution equations, :math:`\psi^{\prime}` is the predicted update without
relaxation, :math:`\psi^{FP}` is the fine data obtained from spatial and temporal interpolation of the
coarse data, and :math:`n` is the minimum number of grid points from a lateral boundary. The specified and
relaxation regions are applied to all dycore variables :math:`\left[\rho \; \rho\Theta \; U\; V\; W \right]`
on the fine mesh.
where :math:`G` is the RHS of the Navier-Stokes equations, :math:`\psi^{*}` is the state variable at the
current RK stage, :math:`\psi^{BDY}` is temporal interpolation of the observational data, :math:`C_{01} = -\ln(0.01) / ({\rm RealWidth - SpecWidth})`
is a constant that ensure the exponential blending function obtains a value of 0.01 at the last relaxation cell,
and :math:`n` is the minimum number of grid points from a lateral boundary.

Sponge zone domain BCs
----------------------
Expand Down
3 changes: 0 additions & 3 deletions Source/ERF_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -534,8 +534,6 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf)

void ERF::init_stuff(int lev, const BoxArray& ba, const DistributionMapping& dm)
{
amrex::Print() << "init_stuff " << ba << std::endl;

if (lev == 0) {
min_k_at_level[lev] = 0;
max_k_at_level[lev] = geom[lev].Domain().bigEnd(2);
Expand All @@ -549,7 +547,6 @@ void ERF::init_stuff(int lev, const BoxArray& ba, const DistributionMapping& dm)
max_k_at_level[lev] = std::max(max_k_at_level[lev], ba[n].bigEnd(2));
}
}
amrex::Print() << "MIN/MAX K AT LEVEL " << lev << " ARE " << min_k_at_level[lev] << " " << max_k_at_level[lev] << std::endl;

// The number of ghost cells for density must be 1 greater than that for velocity
// so that we can go back in forth betwen velocity and momentum on all faces
Expand Down
5 changes: 3 additions & 2 deletions Source/TimeIntegration/ERF_slow_rhs_post.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -472,7 +472,7 @@ void erf_slow_rhs_post (int level, int finest_level,
// cell here if it is present.

// The width to do RHS augmentation
if (width > set_width+1) width -= 1;
if (width > set_width+1) width -= 2;

// Relaxation constants
Real F1 = 1./(10.*dt);
Expand Down Expand Up @@ -571,12 +571,13 @@ void erf_slow_rhs_post (int level, int finest_level,
+ alpha * bdatyhi_np1(ii,jj,k);
});


// NOTE: We pass 'old_cons' here since the tendencies are with
// respect to the start of the RK integration.

// Compute RHS in specified region
//==========================================================
if (set_width > 0 ) {
if (set_width > 0) {
compute_interior_ghost_bxs_xy(tbx, domain, width, 0,
bx_xlo, bx_xhi,
bx_ylo, bx_yhi);
Expand Down
7 changes: 5 additions & 2 deletions Source/Utils/InteriorGhostCells.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -456,13 +456,16 @@ realbdy_compute_interior_ghost_rhs (const std::string& /*init_type*/,
const auto& dom_hi = ubound(domain);
const auto& dom_lo = lbound(domain);

int width2 = width;
if (ivar_idx == IntVars::cons) width2 -= 1;

#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for ( MFIter mfi(S_cur_data[ivar_idx],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
Box tbx = mfi.tilebox();
Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi;
compute_interior_ghost_bxs_xy(tbx, domain, width, set_width,
compute_interior_ghost_bxs_xy(tbx, domain, width2, set_width,
tbx_xlo, tbx_xhi,
tbx_ylo, tbx_yhi);

Expand Down Expand Up @@ -494,7 +497,7 @@ realbdy_compute_interior_ghost_rhs (const std::string& /*init_type*/,
}

wrfbdy_compute_laplacian_relaxation(icomp, 1,
width, set_width, dom_lo, dom_hi, F1, F2,
width2, set_width, dom_lo, dom_hi, F1, F2,
tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi,
arr_xlo, arr_xhi, arr_ylo, arr_yhi,
data_arr, rhs_arr);
Expand Down
22 changes: 13 additions & 9 deletions Source/Utils/Utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,7 +228,7 @@ wrfbdy_compute_laplacian_relaxation (const int& icomp,
int Relax_z = width - Spec_z + 1;
amrex::Real num = amrex::Real(Spec_z + Relax_z);
amrex::Real denom = amrex::Real(Relax_z - 1);
amrex::Real SpecExp = 0.5;
amrex::Real SpecExp = -std::log(0.01) / amrex::Real(width - Spec_z);
amrex::ParallelFor(bx_xlo, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
// Corners with x boxes
Expand All @@ -238,7 +238,7 @@ wrfbdy_compute_laplacian_relaxation (const int& icomp,
int n_ind = std::min(i-dom_lo.x,jj) + 1;
AMREX_ASSERT(n_ind > Spec_z);
amrex::Real Factor = (num - amrex::Real(n_ind))/denom
* std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1));
* std::exp(-SpecExp * amrex::Real(n_ind - Spec_z));
amrex::Real d = data_arr(i ,j ,k ,n+icomp);
amrex::Real d_ip1 = data_arr(i+1,j ,k ,n+icomp);
amrex::Real d_im1 = data_arr(i-1,j ,k ,n+icomp);
Expand All @@ -250,7 +250,8 @@ wrfbdy_compute_laplacian_relaxation (const int& icomp,
amrex::Real delta_yp = arr_xlo(i ,j+1,k,n) - d_jp1;
amrex::Real delta_ym = arr_xlo(i ,j-1,k,n) - d_jm1;
amrex::Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
rhs_arr(i,j,k,n+icomp) += (F1*delta - F2*Laplacian) * Factor;
amrex::Real Temp = (F1*delta - F2*Laplacian) * Factor;
rhs_arr(i,j,k,n+icomp) += Temp;
},
bx_xhi, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
Expand All @@ -261,7 +262,7 @@ wrfbdy_compute_laplacian_relaxation (const int& icomp,
int n_ind = std::min(dom_hi.x-i,jj) + 1;
AMREX_ASSERT(n_ind > Spec_z);
amrex::Real Factor = (num - amrex::Real(n_ind))/denom
* std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1));
* std::exp(-SpecExp * amrex::Real(n_ind - Spec_z));
amrex::Real d = data_arr(i ,j ,k ,n+icomp);
amrex::Real d_ip1 = data_arr(i+1,j ,k ,n+icomp);
amrex::Real d_im1 = data_arr(i-1,j ,k ,n+icomp);
Expand All @@ -273,7 +274,8 @@ wrfbdy_compute_laplacian_relaxation (const int& icomp,
amrex::Real delta_yp = arr_xhi(i ,j+1,k,n) - d_jp1;
amrex::Real delta_ym = arr_xhi(i ,j-1,k,n) - d_jm1;
amrex::Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
rhs_arr(i,j,k,n+icomp) += (F1*delta - F2*Laplacian) * Factor;
amrex::Real Temp = (F1*delta - F2*Laplacian) * Factor;
rhs_arr(i,j,k,n+icomp) += Temp;
});

amrex::ParallelFor(bx_ylo, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -282,7 +284,7 @@ wrfbdy_compute_laplacian_relaxation (const int& icomp,
int n_ind = j - dom_lo.y + 1;
AMREX_ASSERT(n_ind > Spec_z);
amrex::Real Factor = (num - amrex::Real(n_ind))/denom
* std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1));
* std::exp(-SpecExp * amrex::Real(n_ind - Spec_z));
amrex::Real d = data_arr(i ,j ,k ,n+icomp);
amrex::Real d_ip1 = data_arr(i+1,j ,k ,n+icomp);
amrex::Real d_im1 = data_arr(i-1,j ,k ,n+icomp);
Expand All @@ -294,15 +296,16 @@ wrfbdy_compute_laplacian_relaxation (const int& icomp,
amrex::Real delta_yp = arr_ylo(i ,j+1,k,n) - d_jp1;
amrex::Real delta_ym = arr_ylo(i ,j-1,k,n) - d_jm1;
amrex::Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
rhs_arr(i,j,k,n+icomp) += (F1*delta - F2*Laplacian) * Factor;
amrex::Real Temp = (F1*delta - F2*Laplacian) * Factor;
rhs_arr(i,j,k,n+icomp) += Temp;
},
bx_yhi, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
// No corners for y boxes
int n_ind = dom_hi.y - j + 1;
AMREX_ASSERT(n_ind > Spec_z);
amrex::Real Factor = (num - amrex::Real(n_ind))/denom
* std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1));
* std::exp(-SpecExp * amrex::Real(n_ind - Spec_z));
amrex::Real d = data_arr(i ,j ,k ,n+icomp);
amrex::Real d_ip1 = data_arr(i+1,j ,k ,n+icomp);
amrex::Real d_im1 = data_arr(i-1,j ,k ,n+icomp);
Expand All @@ -314,7 +317,8 @@ wrfbdy_compute_laplacian_relaxation (const int& icomp,
amrex::Real delta_yp = arr_yhi(i ,j+1,k,n) - d_jp1;
amrex::Real delta_ym = arr_yhi(i ,j-1,k,n) - d_jm1;
amrex::Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
rhs_arr(i,j,k,n+icomp) += (F1*delta - F2*Laplacian) * Factor;
amrex::Real Temp = (F1*delta - F2*Laplacian) * Factor;
rhs_arr(i,j,k,n+icomp) += Temp;
});
}
#endif

0 comments on commit 2e777a2

Please sign in to comment.