From c402eaeb2e7caa8df4364e648c8572d9336cc5bc Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Sat, 14 Dec 2024 22:42:05 -0500
Subject: [PATCH] update the ambient BCs for compressible this now allows stuff
 to flow out but not in also tweak the initial conditions for convection

---
 pyro/compressible/BC.py                      | 25 ++++++++++++++++----
 pyro/compressible/problems/inputs.convection | 13 +++++-----
 2 files changed, 28 insertions(+), 10 deletions(-)

diff --git a/pyro/compressible/BC.py b/pyro/compressible/BC.py
index 6cdb0eca8..a00efea26 100644
--- a/pyro/compressible/BC.py
+++ b/pyro/compressible/BC.py
@@ -20,7 +20,10 @@
 
 def user(bc_name, bc_edge, variable, ccdata):
     """
-    A hydrostatic boundary.  This integrates the equation of HSE into
+    Extra boundary condition types for compressible hydro.  This includes
+    a hydrostatic boundary, a ramp, and ambient.
+
+    For HSE, this integrates the equation of HSE into
     the ghost cells to get the pressure and density under the assumption
     that the specific internal energy is constant.
 
@@ -152,6 +155,16 @@ def user(bc_name, bc_edge, variable, ccdata):
 
         if bc_edge == "yrb":
 
+            # store the normal momentum -- skip this if we are filling
+            # the sources
+            dens_inside = None
+            mom_normal = None
+            normal_vel_inside = None
+            if "y-momentum" in ccdata.names:
+                mom_normal = ccdata.get_var("y-momentum")[:, myg.jhi]
+                dens_inside = ccdata.get_var("density")[:, myg.jhi]
+                normal_vel_inside = mom_normal / dens_inside
+
             # upper y boundary
 
             # by default, use a zero gradient
@@ -169,13 +182,17 @@ def user(bc_name, bc_edge, variable, ccdata):
 
             elif variable == "y-momentum":
                 rhov = ambient_rho * ambient_v
-                v[:, myg.jhi+1:myg.jhi+myg.ng+1] = rhov
+                # allow stuff to flow out, otherwise, reflect velocity
+                for j in range(myg.jhi+1, myg.jhi+myg.ng+1):
+                    v[:, j] = np.where(mom_normal > 0, mom_normal, rhov)
 
             elif variable == "energy":
                 gamma = ccdata.get_aux("gamma")
-                ke = 0.5 * ambient_rho * (ambient_u**2 + ambient_v**2)
+                yvel = np.where(mom_normal > 0, mom_normal/dens_inside, ambient_v)
+                ke = 0.5 * ambient_rho * (ambient_u**2 + yvel**2)
                 rhoE = ambient_p / (gamma - 1.0) + ke
-                v[:, myg.jhi+1:myg.jhi+myg.ng+1] = rhoE
+                for j in range(myg.jhi+1, myg.jhi+myg.ng+1):
+                    v[:, j] = rhoE[:]
 
         else:
             msg.fail("error: ambient BC not supported for xlb, xrb, or ylb")
diff --git a/pyro/compressible/problems/inputs.convection b/pyro/compressible/problems/inputs.convection
index 150408310..8b337bbf9 100644
--- a/pyro/compressible/problems/inputs.convection
+++ b/pyro/compressible/problems/inputs.convection
@@ -7,17 +7,18 @@ tmax = 10.0
 
 [io]
 basename = convection_
-n_out = 100
+n_out = 100000
+dt_out = 0.25
 
 
 [mesh]
 nx = 128
-ny = 384
+ny = 512
 xmax = 4.0
-ymax = 12.0
+ymax = 16.0
 
-xlboundary = outflow
-xrboundary = outflow
+xlboundary = periodic
+xrboundary = periodic
 
 ylboundary = reflect
 yrboundary = ambient
@@ -26,7 +27,7 @@ yrboundary = ambient
 [convection]
 scale_height = 2.0
 dens_base = 1000.0
-dens_cutoff = 1.e-3
+dens_cutoff = 1.e-4
 
 e_rate = 0.5