fix oops in yesterday's fix (#1957)

erf-model · Nov 20, 2024 · f9ce37e · f9ce37e
1 parent 1a55772
commit f9ce37e
Show file tree

Hide file tree

Showing 3 changed files with 47 additions and 11 deletions.
diff --git a/Source/TimeIntegration/ERF_TI_no_substep_fun.H b/Source/TimeIntegration/ERF_TI_no_substep_fun.H
@@ -12,7 +12,6 @@
         int n_data = IntVars::NumTypes;
 
         const auto& dxInv = fine_geom.InvCellSizeArray();
-        Box domain = fine_geom.Domain();
 
         const amrex::GpuArray<int, IntVars::NumTypes> scomp_fast = {0,0,0,0};
         const amrex::GpuArray<int, IntVars::NumTypes> ncomp_fast = {2,1,1,1};

diff --git a/Source/Utils/ERF_PoissonSolve.cpp b/Source/Utils/ERF_PoissonSolve.cpp
@@ -77,13 +77,30 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
                 }
             });
         } // mfi
-    } // use_terrain
+        //
+        // Note we compute the divergence after we convert rho0W --> Omega
+        //
+        for ( MFIter mfi(rhs[0],TilingIfNotGPU()); mfi.isValid(); ++mfi)
+        {
+            Box bx = mfi.tilebox();
+            const Array4<Real const>& rho0u_arr = mom_mf[IntVars::xmom].const_array(mfi);
+            const Array4<Real const>& rho0v_arr = mom_mf[IntVars::ymom].const_array(mfi);
+            const Array4<Real const>& rho0w_arr = mom_mf[IntVars::zmom].const_array(mfi);
+            const Array4<Real      >&  rhs_arr = rhs[0].array(mfi);
+
+            Real* stretched_dz_d_ptr = stretched_dz_d[lev].data();
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
+                Real dz = stretched_dz_d_ptr[k];
+                rhs_arr(i,j,k) =  (rho0u_arr(i+1,j,k) - rho0u_arr(i,j,k)) * dxInv[0]
+                                 +(rho0v_arr(i,j+1,k) - rho0v_arr(i,j,k)) * dxInv[1]
+                                 +(rho0w_arr(i,j,k+1) - rho0w_arr(i,j,k)) / dz;
+            });
+        } // mfi
 
-    computeDivergence(rhs[0], rho0_u_const, geom_tmp[0]);
+    } else {
+
+        computeDivergence(rhs[0], rho0_u_const, geom_tmp[0]);
 
-    if (l_use_terrain)
-    {
-        MultiFab::Divide(rhs[0],*detJ_cc[0],0,0,1,0);
     }
 
     Real rhsnorm = rhs[0].norm0();
@@ -146,14 +163,34 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     //
     if (mg_verbose > 0)
     {
-        computeDivergence(rhs[0], rho0_u_const, geom_tmp[0]);
         if (l_use_terrain)
         {
-            MultiFab::Divide(rhs[0],*detJ_cc[0],0,0,1,0);
-        } // use_terrain
+            //
+            // Note we compute the divergence before we convert Omega back to rho0W
+            //
+            for ( MFIter mfi(rhs[0],TilingIfNotGPU()); mfi.isValid(); ++mfi)
+            {
+                Box bx = mfi.tilebox();
+                const Array4<Real const>& rho0u_arr = mom_mf[IntVars::xmom].const_array(mfi);
+                const Array4<Real const>& rho0v_arr = mom_mf[IntVars::ymom].const_array(mfi);
+                const Array4<Real const>& rho0w_arr = mom_mf[IntVars::zmom].const_array(mfi);
+                const Array4<Real      >&  rhs_arr = rhs[0].array(mfi);
+
+                Real* stretched_dz_d_ptr = stretched_dz_d[lev].data();
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
+                    Real dz = stretched_dz_d_ptr[k];
+                    rhs_arr(i,j,k) =  (rho0u_arr(i+1,j,k) - rho0u_arr(i,j,k)) * dxInv[0]
+                                     +(rho0v_arr(i,j+1,k) - rho0v_arr(i,j,k)) * dxInv[1]
+                                     +(rho0w_arr(i,j,k+1) - rho0w_arr(i,j,k)) / dz;
+                });
+            } // mfi
+
+        } else {
+            computeDivergence(rhs[0], rho0_u_const, geom_tmp[0]);
+        }
 
         amrex::Print() << "Max norm of divergence after solve at level " << lev << " : " << rhs[0].norm0() << std::endl;
-    }
+    } // mg_verbose
 
     //
     // ****************************************************************************

diff --git a/Source/Utils/ERF_solve_with_fft.cpp b/Source/Utils/ERF_solve_with_fft.cpp
@@ -233,7 +233,7 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array<MultiFab,
                     fz_arr(i,j,k) = -(p_arr(i,j,k) - p_arr(i,j,k-1)) / dz;
                 }
             });
-        } else {
+        } else { // no grid stretching
             const Real dz_inv = dxInv[2];
             ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {