make w2a initialization more accurate, especially for low resolutions (#1669)

* make w2a initialization more accurate, especially for low resolutions

* formatting

* avoid dividing by 0, just in case

Author: mbkuhn

amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H

@@ -652,6 +652,16 @@ struct InitDataOp<W2AWaves>
         const bool has_beach = wdata.has_beach && multiphase_mode;
         const bool init_wave_field = wdata.init_wave_field || !multiphase_mode;
 
+        amrex::MultiFab phi_tmp_fab(
+            velocity(level).boxArray(), velocity(level).DistributionMap(), 1,
+            3);
+        amrex::MultiFab vel_tmp_fab(
+            velocity(level).boxArray(), velocity(level).DistributionMap(),
+            AMREX_SPACEDIM, 3);
+
+        const auto& phi_tmp = phi_tmp_fab.arrays();
+        const auto& vel_tmp = vel_tmp_fab.arrays();
+
         amrex::ParallelFor(
             velocity(level), amrex::IntVect(3),
             [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
@@ -673,17 +683,13 @@ struct InitDataOp<W2AWaves>
                     x, problo[0], gen_length, probhi[0], beach_length, wave_sol,
                     bulk, outlet);
 
-                const amrex::Real phi = local_profile[3] - z;
-                const amrex::Real cell_length_2D =
-                    std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]);
-                if (phi + cell_length_2D >= 0) {
-                    vel[nbx](i, j, k, 0) = local_profile[0];
-                    vel[nbx](i, j, k, 1) = local_profile[1];
-                    vel[nbx](i, j, k, 2) = local_profile[2];
-                }
+                phi_tmp[nbx](i, j, k) = local_profile[3] - z;
+                vel_tmp[nbx](i, j, k, 0) = local_profile[0];
+                vel_tmp[nbx](i, j, k, 1) = local_profile[1];
+                vel_tmp[nbx](i, j, k, 2) = local_profile[2];
 
                 if (multiphase_mode) {
-                    phi_arrs[nbx](i, j, k) = phi;
+                    phi_arrs[nbx](i, j, k) = phi_tmp[nbx](i, j, k);
                 }
 
                 // Default w2a values matter for where no updates happen
@@ -693,6 +699,39 @@ struct InitDataOp<W2AWaves>
                 w2a_vel[nbx](i, j, k, 2) = quiescent[0];
             });
         amrex::Gpu::streamSynchronize();
+
+        amrex::ParallelFor(
+            velocity(level), amrex::IntVect(2),
+            [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+                const amrex::Real eps = 2. * std::cbrt(dx[0] * dx[1] * dx[2]);
+                const amrex::Real vof_local =
+                    multiphase::levelset_to_vof(i, j, k, eps, phi_tmp[nbx]);
+
+                if (vof_local > 1. - constants::TIGHT_TOL) {
+                    // Entire cell is liquid
+                    vel[nbx](i, j, k, 0) = vel_tmp[nbx](i, j, k, 0);
+                    vel[nbx](i, j, k, 1) = vel_tmp[nbx](i, j, k, 1);
+                    vel[nbx](i, j, k, 2) = vel_tmp[nbx](i, j, k, 2);
+                } else if (vof_local > 1e-3) {
+                    // Velocity of liquid becomes velocity of entire cell
+                    // Interpolate using cell-centered values to liquid centroid
+                    // Consider only z-direction to find centroid
+                    const amrex::Real z_cell = problo[2] + (k + 0.5) * dx[2];
+                    const amrex::Real z_cell_below = z_cell - dx[2];
+                    const amrex::Real z_cell_bottom = z_cell - 0.5 * dx[2];
+                    const amrex::Real z_liq_center =
+                        z_cell_bottom + 0.5 * vof_local * dx[2];
+                    for (int n = 0; n < AMREX_SPACEDIM; ++n) {
+                        vel[nbx](i, j, k, n) =
+                            vel_tmp[nbx](i, j, k - 1, n) +
+                            (vel_tmp[nbx](i, j, k, n) -
+                             vel_tmp[nbx](i, j, k - 1, n)) *
+                                ((z_liq_center - z_cell_below) /
+                                 (z_cell - z_cell_below + constants::EPS));
+                    }
+                }
+            });
+        amrex::Gpu::streamSynchronize();
 #else
         amrex::ignore_unused(data, level, geom, multiphase_mode);
 #endif