From 2e777a25eb524663dcc1e7eb009645566fe0ab0d Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Wed, 13 Mar 2024 15:19:35 -0700 Subject: [PATCH] Nudging Mods (#1492) * Best nudging FLOWMAS simulation yet. * Clean up and add doc description. * Clarify docs. --------- Co-authored-by: Aaron Lattanzi Co-authored-by: Ann Almgren --- Docs/sphinx_doc/BoundaryConditions.rst | 34 +++++++++++--------- Source/ERF_make_new_level.cpp | 3 -- Source/TimeIntegration/ERF_slow_rhs_post.cpp | 5 +-- Source/Utils/InteriorGhostCells.cpp | 7 ++-- Source/Utils/Utils.H | 22 +++++++------ 5 files changed, 39 insertions(+), 32 deletions(-) diff --git a/Docs/sphinx_doc/BoundaryConditions.rst b/Docs/sphinx_doc/BoundaryConditions.rst index 8e835b8cd..1d3c81864 100644 --- a/Docs/sphinx_doc/BoundaryConditions.rst +++ b/Docs/sphinx_doc/BoundaryConditions.rst @@ -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 = `` (yellow + blue) and analogously the width of the interior Dirichlet region may be specified with ``erf.real_set_width = `` (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 @@ -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 ---------------------- diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index c23c6199b..092c2c898 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -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); @@ -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 diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index 187fb466c..7e2b130da 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -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); @@ -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); diff --git a/Source/Utils/InteriorGhostCells.cpp b/Source/Utils/InteriorGhostCells.cpp index 09872574d..07e201962 100644 --- a/Source/Utils/InteriorGhostCells.cpp +++ b/Source/Utils/InteriorGhostCells.cpp @@ -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); @@ -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); diff --git a/Source/Utils/Utils.H b/Source/Utils/Utils.H index ddc6607bf..0058f085f 100644 --- a/Source/Utils/Utils.H +++ b/Source/Utils/Utils.H @@ -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 @@ -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); @@ -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 { @@ -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); @@ -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 @@ -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); @@ -294,7 +296,8 @@ 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 { @@ -302,7 +305,7 @@ wrfbdy_compute_laplacian_relaxation (const int& icomp, 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); @@ -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