Fixes for 2d squall line simulation - Introducing new bc - HO_Outflow (…

…#1602) * Hard code lin extrap for scalars at top * Generalizing linextrap * Generalizing linextrap bc * Generalizing linextrap final * Finalizing generalization of linextrap * Checking with some changes * Checking with some changes * Box issues in ERFPhysBCFunct_cons::operator() fixed * hoextrapcc for all variables * Removing unwanted comments * New bc HO_outflow * fast coefficient change for ho_outflow as well * Correct link for SAM * Updating inputs file with zhi bc as HO_Outflow --------- Co-authored-by: Mahesh Natarajan <[email protected]> Co-authored-by: Aaron M. Lattanzi <[email protected]>
erf-model · May 6, 2024 · 8f721f5 · 8f721f5
1 parent cbb3fff
commit 8f721f5
Show file tree

Hide file tree

Showing 9 changed files with 54 additions and 12 deletions.
diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst
@@ -119,7 +119,7 @@ The total production rates including the contribution from aggregation, accretio
 bergeron process, freezing and autoconversion are listed below without derivation.
 For details, please refer to Yuh-Lang Lin et al (J. Climate Appl. Meteor, 22, 1065, 1983) and
 Marat F. Khairoutdinov and David A. Randall's (J. Atm Sciences, 607, 1983).
-The implementation of microphysics model in ERF is similar to the that in the SAM code (http://rossby.msrc.sunysb.edu/~marat/SAM.html)
+The implementation of microphysics model in ERF is similar to the that in the SAM code (http://rossby.msrc.sunysb.edu/SAM.html)
 
 Accretion
 ~~~~~~~~~~~~

diff --git a/Exec/SquallLine_2D/inputs_moisture_Gabersek b/Exec/SquallLine_2D/inputs_moisture_Gabersek
@@ -15,9 +15,8 @@ geometry.is_periodic = 0 1 0
 
 xlo.type = "Open"
 xhi.type = "Open"
-
 zlo.type = "SlipWall"
-zhi.type = "Outflow"
+zhi.type = "HO_Outflow"
 
 # TIME STEP CONTROL
 erf.use_native_mri = 1

diff --git a/Exec/SquallLine_2D/inputs_moisture_WRF b/Exec/SquallLine_2D/inputs_moisture_WRF
@@ -16,7 +16,7 @@ geometry.is_periodic = 0 1 0
 xlo.type = "Open"
 xhi.type = "Open"
 zlo.type = "SlipWall"
-zhi.type = "Outflow"
+zhi.type = "HO_Outflow"
 
 # TIME STEP CONTROL
 erf.use_native_mri = 1

diff --git a/Source/BoundaryConditions/BoundaryConditions_cons.cpp b/Source/BoundaryConditions/BoundaryConditions_cons.cpp
@@ -108,6 +108,8 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4<Real>& dest_arr,
                     dest_arr(i,j,k,icomp+n) =  dest_arr(iflip,j,k,icomp+n);
                 } else if (bc_ptr[icomp+n].lo(0) == ERFBCType::reflect_odd) {
                     dest_arr(i,j,k,icomp+n) = -dest_arr(iflip,j,k,icomp+n);
+                } else if (bc_ptr[icomp+n].lo(0) == ERFBCType::hoextrapcc) {
+                    dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(dom_lo.x,j,k,icomp+n) - dest_arr(dom_lo.x+1,j,k,icomp+n) ;
                 }
             },
             bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
@@ -120,6 +122,8 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4<Real>& dest_arr,
                     dest_arr(i,j,k,icomp+n) =  dest_arr(iflip,j,k,icomp+n);
                 } else if (bc_ptr[icomp+n].hi(0) == ERFBCType::reflect_odd) {
                     dest_arr(i,j,k,icomp+n) = -dest_arr(iflip,j,k,icomp+n);
+                } else if (bc_ptr[icomp+n].hi(0) == ERFBCType::hoextrapcc) {
+                    dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(dom_hi.x,j,k,icomp+n) - dest_arr(dom_hi.x-1,j,k,icomp+n) ;
                 }
             }
         );
@@ -141,7 +145,10 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4<Real>& dest_arr,
                     dest_arr(i,j,k,icomp+n) =  dest_arr(i,jflip,k,icomp+n);
                 } else if (bc_ptr[icomp+n].lo(1) == ERFBCType::reflect_odd) {
                     dest_arr(i,j,k,icomp+n) = -dest_arr(i,jflip,k,icomp+n);
+                } else if (bc_ptr[icomp+n].lo(1) == ERFBCType::hoextrapcc) {
+                    dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(i,dom_lo.y,k,icomp+n) - dest_arr(i,dom_lo.y+1,k,icomp+n) ;
                 }
+
             },
             bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
                 int jflip =  2*dom_hi.y + 1 - j;
@@ -153,6 +160,8 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4<Real>& dest_arr,
                     dest_arr(i,j,k,icomp+n) =  dest_arr(i,jflip,k,icomp+n);
                 } else if (bc_ptr[icomp+n].hi(1) == ERFBCType::reflect_odd) {
                     dest_arr(i,j,k,icomp+n) = -dest_arr(i,jflip,k,icomp+n);
+                } else if (bc_ptr[icomp+n].hi(1) == ERFBCType::hoextrapcc) {
+                    dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(i,dom_hi.y,k,icomp+n) - dest_arr(i,dom_hi.y-1,k,icomp+n);
                 }
             }
         );
@@ -249,6 +258,8 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& dest_arr
                     Real delta_z = (dom_lo.z - k) / dxInv[2];
                     dest_arr(i,j,k,icomp+n) = dest_arr(i,j,dom_lo.z,icomp+n) -
                         delta_z*l_bc_neumann_vals_d[icomp+n][2]*dest_arr(i,j,dom_lo.z,Rho_comp);
+                } else if (bc_ptr[icomp+n].lo(2) == ERFBCType::hoextrapcc) {
+                    dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(i,j,dom_lo.z,icomp+n) - dest_arr(i,j,dom_lo.z+1,icomp+n);
                 }
             },
             bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
@@ -270,6 +281,8 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& dest_arr
                         dest_arr(i,j,k,icomp+n) = dest_arr(i,j,dom_hi.z,icomp+n) +
                             delta_z*l_bc_neumann_vals_d[icomp+n][5]*dest_arr(i,j,dom_hi.z,Rho_comp);
                     }
+                } else if (bc_ptr[icomp+n].hi(2) == ERFBCType::hoextrapcc){
+                    dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(i,j,dom_hi.z,icomp+n) - dest_arr(i,j,dom_hi.z-1,icomp+n);
                 }
             }
         );

diff --git a/Source/IndexDefines.H b/Source/IndexDefines.H
@@ -136,7 +136,7 @@ namespace EddyDiff {
 }
 
 enum struct ERF_BC {
-    symmetry, inflow, outflow, open, no_slip_wall, slip_wall, periodic, MOST, undefined
+    symmetry, inflow, outflow, ho_outflow, open, no_slip_wall, slip_wall, periodic, MOST, undefined
 };
 
 // NOTE: the first of these must match up with the BCType enum
@@ -154,7 +154,8 @@ enum mathematicalBndryTypes : int {
     ext_dir_ingested =  102,
     neumann          =  103,
     neumann_int      =  104,
-    open             =  105
+    open             =  105,
+    hoextrapcc         =  106
 };
 }
 

diff --git a/Source/Initialization/ERF_init_bcs.cpp b/Source/Initialization/ERF_init_bcs.cpp
@@ -80,6 +80,12 @@ void ERF::init_bcs ()
               phys_bc_type[ori] = ERF_BC::open;
             domain_bc_type[ori] = "Open";
         }
+        else if (bc_type == "ho_outflow")
+        {
+            phys_bc_type[ori] = ERF_BC::ho_outflow;
+            domain_bc_type[ori] = "HO_Outflow";
+        }
+
         else if (bc_type == "inflow")
         {
             // Print() << bcid << " set to inflow.\n";
@@ -283,7 +289,7 @@ void ERF::init_bcs ()
                     domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::reflect_odd);
                 }
             }
-            else if (bct == ERF_BC::outflow)
+            else if (bct == ERF_BC::outflow or bct == ERF_BC::ho_outflow )
             {
                 if (side == Orientation::low) {
                     for (int i = 0; i < AMREX_SPACEDIM; i++) {
@@ -412,6 +418,18 @@ void ERF::init_bcs ()
                     }
                 }
             }
+            else if ( bct == ERF_BC::ho_outflow )
+            {
+                if (side == Orientation::low) {
+                    for (int i = 0; i < NVAR_max; i++) {
+                        domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::hoextrapcc);
+                    }
+                } else {
+                    for (int i = 0; i < NVAR_max; i++) {
+                        domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::hoextrapcc);
+                    }
+                }
+            }
             else if ( bct == ERF_BC::open )
             {
                 if (side == Orientation::low) {

diff --git a/Source/TimeIntegration/ERF_make_fast_coeffs.cpp b/Source/TimeIntegration/ERF_make_fast_coeffs.cpp
@@ -35,7 +35,8 @@ void make_fast_coeffs (int /*level*/,
                        Real gravity, Real c_p,
                        std::unique_ptr<MultiFab>& detJ_cc,
                        const MultiFab* r0, const MultiFab* pi0,
-                       Real dtau, Real beta_s)
+                       Real dtau, Real beta_s,
+                       amrex::GpuArray<ERF_BC, AMREX_SPACEDIM*2> &phys_bc_type)
 {
     BL_PROFILE_VAR("make_fast_coeffs()",make_fast_coeffs);
 
@@ -215,6 +216,10 @@ void make_fast_coeffs (int /*level*/,
 
           // w_khi = 0
           coeffA_a(i,j,hi.z+1) =  0.0;
+          if(phys_bc_type[5] == ERF_BC::outflow or phys_bc_type[5] == ERF_BC::ho_outflow)
+          {
+              coeffA_a(i,j,hi.z+1) =  -1.0;
+          }
           coeffB_a(i,j,hi.z+1) =  1.0;
           coeffC_a(i,j,hi.z+1) =  0.0;
 
@@ -239,7 +244,12 @@ void make_fast_coeffs (int /*level*/,
         for (int j = lo.y; j <= hi.y; ++j) {
             AMREX_PRAGMA_SIMD
             for (int i = lo.x; i <= hi.x; ++i) {
+
                 coeffA_a(i,j,hi.z+1) =  0.0;
+                  if(phys_bc_type[5] == ERF_BC::outflow or phys_bc_type[5] == ERF_BC::ho_outflow)
+                {
+                    coeffA_a(i,j,hi.z+1) =  -1.0;
+                }
                 coeffB_a(i,j,hi.z+1) =  1.0;
                 coeffC_a(i,j,hi.z+1) =  0.0;
             }

diff --git a/Source/TimeIntegration/TI_fast_headers.H b/Source/TimeIntegration/TI_fast_headers.H
@@ -118,6 +118,7 @@ void make_fast_coeffs (int level,
                        const amrex::MultiFab* r0,
                        const amrex::MultiFab* pi0,
                        const amrex::Real dtau,
-                       const amrex::Real beta_s);
+                       const amrex::Real beta_s,
+                       amrex::GpuArray<ERF_BC, AMREX_SPACEDIM*2> &phys_bc_type);
 
 #endif
diff --git a/Source/TimeIntegration/TI_fast_rhs_fun.H b/Source/TimeIntegration/TI_fast_rhs_fun.H
@@ -88,7 +88,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
             // Note we pass in the *old* detJ here
             make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom,
                              l_use_moisture, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p,
-                             detJ_cc[level], r0, pi0, dtau, beta_s);
+                             detJ_cc[level], r0, pi0, dtau, beta_s, phys_bc_type);
 
             if (fast_step == 0) {
                 // If this is the first substep we pass in S_old as the previous step's solution
@@ -121,7 +121,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
                 // If this is the first substep we make the coefficients since they are based only on stage data
                 make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom,
                                  l_use_moisture, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p,
-                                 detJ_cc[level], r0, pi0, dtau, beta_s);
+                                 detJ_cc[level], r0, pi0, dtau, beta_s, phys_bc_type);
 
                 // If this is the first substep we pass in S_old as the previous step's solution
                 erf_fast_rhs_T(fast_step, nrk, level, finest_level,
@@ -145,7 +145,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
                 // If this is the first substep we make the coefficients since they are based only on stage data
                 make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom,
                                  l_use_moisture, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p,
-                                 detJ_cc[level], r0, pi0, dtau, beta_s);
+                                 detJ_cc[level], r0, pi0, dtau, beta_s, phys_bc_type);
 
                 // If this is the first substep we pass in S_old as the previous step's solution
                 erf_fast_rhs_N(fast_step, nrk, level, finest_level,