Wave energy postprocessing routine (Exawind#615)

* setting up postprocessing routine to output wave energy, both kinetic and potential. Assumes that gravity acts in the negative z direction.

* format and cleanup

* changes for CI

* more modifications for CI

* Updating the wave energy function to do better in areas where there is an air bubble between some liquid and the main surface. Also, comparing to the exact solution in the unit test.

amr-wind/utilities/sampling/CMakeLists.txt

@@ -12,4 +12,5 @@ target_sources(${amr_wind_lib_name}
   KineticEnergy.cpp
   Enstrophy.cpp
   FreeSurface.cpp
+  WaveEnergy.cpp
   )

amr-wind/utilities/sampling/WaveEnergy.H

@@ -0,0 +1,106 @@
+#ifndef WAVEENERGY_H
+#define WAVEENERGY_H
+
+#include "amr-wind/CFDSim.H"
+#include "amr-wind/utilities/PostProcessing.H"
+
+namespace amr_wind {
+namespace wave_energy {
+
+/** wave energy object
+ *  \ingroup Wave Energy
+ *
+ *  A concrete implementation of the post-processing interface that deals with
+ *  integrating total wave energy. This routine calculates the kinetic energy
+ *  (KE) and potential energy (PE) according to the following definitions. It
+ *  allows for a user-defined potential energy offset, as well (PE_0), intended
+ *  to give the result PE = 0 for an undisturbed interface.
+ *
+ *  KE = 0.5 * (Integral over liquid phase) \rho * u^2 dV
+ *  PE = (Integral over liquid phase) \rho * g * z_liq dV + PE_0
+ *
+ *  Note: the output quantities are not normalized.
+ */
+class WaveEnergy : public PostProcessBase::Register<WaveEnergy>
+{
+public:
+    static std::string identifier() { return "WaveEnergy"; }
+
+    WaveEnergy(CFDSim& /*sim*/, std::string /*label*/);
+
+    ~WaveEnergy() override;
+
+    //! Perform actions before mesh is created
+    void pre_init_actions() override {}
+
+    //! Read user inputs and get information needed for calculations
+    void initialize() override;
+
+    //! Integrate energy components and output to file
+    void post_advance_work() override;
+
+    void post_regrid_actions() override {}
+
+    //! Calculate the sum of stated energy in liquid phase
+    amrex::Real calculate_kinetic_energy();
+    amrex::Real calculate_potential_energy();
+
+    //! Output private variables that store energy measurements
+    void wave_energy(amrex::Real& ke, amrex::Real& pe) const
+    {
+        ke = m_wave_kinetic_energy;
+        pe = m_wave_potential_energy;
+    };
+
+private:
+    //! prepare ASCII file and directory
+    void prepare_ascii_file();
+
+    //! Output sampled data in ASCII format
+    virtual void write_ascii();
+
+    //! store the total wave energy
+    amrex::Real m_wave_kinetic_energy{0.0};
+    amrex::Real m_wave_potential_energy{0.0};
+
+    //! Reference to the CFD sim
+    CFDSim& m_sim;
+
+    /** Name of this sampling object.
+     *
+     *  The label is used to read user inputs from file and is also used for
+     *  naming files directories depending on the output format.
+     */
+    const std::string m_label;
+
+    //! reference to velocity
+    const Field& m_velocity;
+    //! reference to vof
+    const Field& m_vof;
+
+    //! gravity vector
+    amrex::Vector<amrex::Real> m_gravity{{0.0, 0.0, -9.81}};
+
+    //! offset for potential energy calculation
+    amrex::Real m_pe_off = 0.0;
+
+    //! density of liquid phase
+    amrex::Real m_rho1 = 10.0;
+
+    //! filename for ASCII output
+    std::string m_out_fname;
+
+    //! Frequency of data sampling and output
+    int m_out_freq{10};
+
+    //! width in ASCII output
+    int m_width{22};
+
+    //! precision in ASCII output
+    int m_precision{12};
+};
+
+} // namespace wave_energy
+} // namespace amr_wind
+
+#endif /* WAVEENERGY_H */

amr-wind/utilities/sampling/WaveEnergy.cpp

@@ -0,0 +1,217 @@
+#include "amr-wind/utilities/sampling/WaveEnergy.H"
+#include "amr-wind/utilities/io_utils.H"
+#include "amr-wind/utilities/ncutils/nc_interface.H"
+#include "amr-wind/physics/multiphase/MultiPhase.H"
+#include <AMReX_MultiFabUtil.H>
+#include <utility>
+#include "AMReX_ParmParse.H"
+#include "amr-wind/utilities/IOManager.H"
+
+namespace amr_wind {
+namespace wave_energy {
+
+WaveEnergy::WaveEnergy(CFDSim& sim, std::string label)
+    : m_sim(sim)
+    , m_label(std::move(label))
+    , m_velocity(sim.repo().get_field("velocity"))
+    , m_vof(sim.repo().get_field("vof"))
+{}
+
+WaveEnergy::~WaveEnergy() = default;
+
+void WaveEnergy::initialize()
+{
+    BL_PROFILE("amr-wind::WaveEnergy::initialize");
+    {
+        amrex::ParmParse pp(m_label);
+        pp.query("output_frequency", m_out_freq);
+        pp.query("potential_energy_offset", m_pe_off);
+    }
+    // Get gravity and density
+    {
+        amrex::ParmParse pp("incflo");
+        pp.queryarr("gravity", m_gravity);
+    }
+    auto& mphase = m_sim.physics_manager().get<amr_wind::MultiPhase>();
+    m_rho1 = mphase.rho1();
+
+    prepare_ascii_file();
+}
+
+amrex::Real WaveEnergy::calculate_kinetic_energy()
+{
+    BL_PROFILE("amr-wind::WaveEnergy::calculate_kinetic_energy");
+
+    // integrated total wave Energy
+    amrex::Real wave_ke = 0.0;
+
+    // Liquid density
+    const amrex::Real rho_liq = m_rho1;
+
+    const int finest_level = m_velocity.repo().num_active_levels() - 1;
+    const auto& geom = m_velocity.repo().mesh().Geom();
+
+    for (int lev = 0; lev <= finest_level; lev++) {
+
+        amrex::iMultiFab level_mask;
+        if (lev < finest_level) {
+            level_mask = makeFineMask(
+                m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
+                m_sim.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
+        } else {
+            level_mask.define(
+                m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
+                1, 0, amrex::MFInfo());
+            level_mask.setVal(1);
+        }
+
+        const amrex::Real cell_vol = geom[lev].CellSize()[0] *
+                                     geom[lev].CellSize()[1] *
+                                     geom[lev].CellSize()[2];
+
+        wave_ke += amrex::ReduceSum(
+            m_vof(lev), m_velocity(lev), level_mask, 0,
+            [=] AMREX_GPU_HOST_DEVICE(
+                amrex::Box const& bx,
+                amrex::Array4<amrex::Real const> const& vof_arr,
+                amrex::Array4<amrex::Real const> const& vel_arr,
+                amrex::Array4<int const> const& mask_arr) -> amrex::Real {
+                amrex::Real Wave_Energy_Fab = 0.0;
+
+                amrex::Loop(
+                    bx, [=, &Wave_Energy_Fab](int i, int j, int k) noexcept {
+                        Wave_Energy_Fab +=
+                            cell_vol * mask_arr(i, j, k) * 0.5 * rho_liq *
+                            vof_arr(i, j, k) *
+                            (vel_arr(i, j, k, 0) * vel_arr(i, j, k, 0) +
+                             vel_arr(i, j, k, 1) * vel_arr(i, j, k, 1) +
+                             vel_arr(i, j, k, 2) * vel_arr(i, j, k, 2));
+                    });
+                return Wave_Energy_Fab;
+            });
+    }
+
+    amrex::ParallelDescriptor::ReduceRealSum(wave_ke);
+
+    return wave_ke;
+}
+
+amrex::Real WaveEnergy::calculate_potential_energy()
+{
+    BL_PROFILE("amr-wind::WaveEnergy::calculate_potential_energy");
+
+    // integrated total wave Energy
+    amrex::Real wave_pe = 0.0;
+    // Liquid density
+    const amrex::Real rho_liq = m_rho1;
+    // gravity constant
+    const amrex::Real g = -m_gravity[2];
+
+    const int finest_level = m_velocity.repo().num_active_levels() - 1;
+    const auto& geom = m_velocity.repo().mesh().Geom();
+
+    for (int lev = 0; lev <= finest_level; lev++) {
+
+        amrex::iMultiFab level_mask;
+        if (lev < finest_level) {
+            level_mask = makeFineMask(
+                m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
+                m_sim.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
+        } else {
+            level_mask.define(
+                m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
+                1, 0, amrex::MFInfo());
+            level_mask.setVal(1);
+        }
+
+        const amrex::Real cell_vol = geom[lev].CellSize()[0] *
+                                     geom[lev].CellSize()[1] *
+                                     geom[lev].CellSize()[2];
+
+        const amrex::Real dz = geom[lev].CellSize()[2];
+        const amrex::Real probloz = geom[lev].ProbLo()[2];
+
+        wave_pe += amrex::ReduceSum(
+            m_vof(lev), level_mask, 0,
+            [=] AMREX_GPU_HOST_DEVICE(
+                amrex::Box const& bx,
+                amrex::Array4<amrex::Real const> const& vof_arr,
+                amrex::Array4<int const> const& mask_arr) -> amrex::Real {
+                amrex::Real Wave_Energy_Fab = 0.0;
+
+                amrex::Loop(
+                    bx, [=, &Wave_Energy_Fab](int i, int j, int k) noexcept {
+                        // Crude model of liquid height in multiphase cells
+                        amrex::Real kk =
+                            (vof_arr(i, j, k + 1) > vof_arr(i, j, k)) ? k + 1
+                                                                      : k;
+                        amrex::Real dir =
+                            (vof_arr(i, j, k + 1) > vof_arr(i, j, k)) ? -1 : 1;
+                        const amrex::Real zl =
+                            probloz + (kk + dir * 0.5 * vof_arr(i, j, k)) * dz;
+                        Wave_Energy_Fab += cell_vol * mask_arr(i, j, k) *
+                                           rho_liq * vof_arr(i, j, k) * g * zl;
+                    });
+                return Wave_Energy_Fab;
+            });
+    }
+
+    amrex::ParallelDescriptor::ReduceRealSum(wave_pe);
+
+    return wave_pe;
+}
+
+void WaveEnergy::post_advance_work()
+{
+    BL_PROFILE("amr-wind::WaveEnergy::post_advance_work");
+    const auto& time = m_sim.time();
+    const int tidx = time.time_index();
+    // Skip processing if it is not an output timestep
+    if (!(tidx % m_out_freq == 0)) {
+        return;
+    }
+
+    m_wave_kinetic_energy = calculate_kinetic_energy();
+    m_wave_potential_energy = calculate_potential_energy() + m_pe_off;
+
+    write_ascii();
+}
+
+void WaveEnergy::prepare_ascii_file()
+{
+    BL_PROFILE("amr-wind::WaveEnergy::prepare_ascii_file");
+
+    const std::string post_dir = "post_processing";
+    const std::string sname =
+        amrex::Concatenate(m_label, m_sim.time().time_index());
+
+    if (!amrex::UtilCreateDirectory(post_dir, 0755)) {
+        amrex::CreateDirectoryFailed(post_dir);
+    }
+    m_out_fname = post_dir + "/" + sname + ".txt";
+
+    if (amrex::ParallelDescriptor::IOProcessor()) {
+        std::ofstream f(m_out_fname.c_str());
+        f << "time_step time kinetic_energy potential_energy" << std::endl;
+        f.close();
+    }
+}
+
+void WaveEnergy::write_ascii()
+{
+    BL_PROFILE("amr-wind::WaveEnergy::write_ascii");
+
+    if (amrex::ParallelDescriptor::IOProcessor()) {
+        std::ofstream f(m_out_fname.c_str(), std::ios_base::app);
+        f << m_sim.time().time_index() << std::scientific
+          << std::setprecision(m_precision) << std::setw(m_width)
+          << m_sim.time().new_time();
+        f << std::setw(m_width) << m_wave_kinetic_energy;
+        f << std::setw(m_width) << m_wave_potential_energy;
+        f << std::endl;
+        f.close();
+    }
+}
+
+} // namespace wave_energy
+} // namespace amr_wind

unit_tests/utilities/CMakeLists.txt

@@ -7,6 +7,7 @@ target_sources(${amr_wind_unit_test_exe_name} PRIVATE
   test_sampling.cpp
   test_linear_interpolation.cpp
   test_free_surface.cpp
+  test_wave_energy.cpp
   )
 
 if (AMR_WIND_ENABLE_NETCDF)