Reduce memory usage of native sampling (Exawind#1207)

.gitignore

@@ -22,3 +22,4 @@ settings.json
 .ccls
 .ccls-cache/
 build*
+__pycache__/

amr-wind/utilities/sampling/Sampling.H

@@ -74,10 +74,15 @@ public:
         return static_cast<int>(m_total_particles);
     }
 
-    int num_output_particles() const
+#ifdef AMR_WIND_USE_NETCDF
+    int num_netcdf_output_particles() const
     {
-        return static_cast<int>(m_output_particles);
+        if (m_out_fmt == "netcdf") {
+            return static_cast<int>(m_netcdf_output_particles);
+        }
+        return 0;
     }
+#endif
 
     const amrex::Vector<std::string>& var_names() const { return m_var_names; }
 
@@ -128,16 +133,19 @@ private:
      */
     const std::string m_label;
 
+    //! Format of the data output (native, ascii, netcdf, etc.)
+#ifdef AMR_WIND_USE_NETCDF
+    std::string m_out_fmt{"netcdf"};
+    std::string m_ncfile_name;
+
     // Current sample buffer for full container
     std::vector<double> m_sample_buf;
 
     // Current sample buffer for output
     std::vector<double> m_output_buf;
 
-    //! Format of the data output (native, ascii, netcdf, etc.)
-#ifdef AMR_WIND_USE_NETCDF
-    std::string m_out_fmt{"netcdf"};
-    std::string m_ncfile_name;
+    //! Number of output particles in netcdf
+    size_t m_netcdf_output_particles{0};
 #else
     std::string m_out_fmt{"native"};
 #endif
@@ -154,9 +162,6 @@ private:
     //! Number of sampling container particles:
     size_t m_total_particles{0};
 
-    //! Number of output particles
-    size_t m_output_particles{0};
-
     //! Frequency of data sampling and output
     int m_out_freq{100};
 

amr-wind/utilities/sampling/Sampling.cpp

@@ -94,16 +94,24 @@ void Sampling::initialize()
 
     update_container();
 
+#ifdef AMR_WIND_USE_NETCDF
     if (m_out_fmt == "netcdf") {
         prepare_netcdf_file();
+        m_sample_buf.assign(m_total_particles * m_var_names.size(), 0.0);
     }
-
-    m_sample_buf.assign(m_total_particles * m_var_names.size(), 0.0);
+#endif
 
     if (m_restart_sample) {
         sampling_workflow();
         sampling_post();
     }
+
+    // Check
+    for (const auto& obj : m_samplers) {
+        if ((obj->do_convert_velocity_los()) && (m_out_fmt != "netcdf")) {
+            amrex::Abort("Velocity line of sight capability requires NetCDF");
+        }
+    }
 }
 
 void Sampling::update_container()
@@ -187,7 +195,11 @@ void Sampling::sampling_post()
         obj->post_sample_actions();
     }
 
-    m_output_buf.clear();
+#ifdef AMR_WIND_USE_NETCDF
+    if (m_out_fmt == "netcdf") {
+        m_output_buf.clear();
+    }
+#endif
 }
 
 void Sampling::post_regrid_actions()
@@ -205,6 +217,11 @@ void Sampling::convert_velocity_lineofsight()
 {
     BL_PROFILE("amr-wind::Sampling::convert_velocity_lineofsight");
 
+    if (m_out_fmt != "netcdf") {
+        return;
+    }
+
+#ifdef AMR_WIND_USE_NETCDF
     amrex::Vector<int> vel_map(AMREX_SPACEDIM, 0);
     const amrex::Vector<std::string> vnames = {
         "velocityx", "velocityy", "velocityz"};
@@ -251,11 +268,22 @@ void Sampling::convert_velocity_lineofsight()
         }
         soffset += sample_size;
     }
+#else
+    amrex::Abort(
+        "NetCDF support was not enabled during build time. Please recompile or "
+        "use native format");
+#endif
 }
 
 void Sampling::create_output_buffer()
 {
     BL_PROFILE("amr-wind::Sampling::create_output_buffer");
+
+    if (m_out_fmt != "netcdf") {
+        return;
+    }
+
+#ifdef AMR_WIND_USE_NETCDF
     const long nvars = m_var_names.size();
     for (int iv = 0; iv < nvars; ++iv) {
         long offset = iv * m_scontainer->num_sampling_particles();
@@ -281,13 +309,27 @@ void Sampling::create_output_buffer()
         }
     }
 
-    m_output_particles = m_output_buf.size() / nvars;
+    m_netcdf_output_particles = m_output_buf.size() / nvars;
+#else
+    amrex::Abort(
+        "NetCDF support was not enabled during build time. Please recompile or "
+        "use native format");
+#endif
 }
 
 void Sampling::fill_buffer()
 {
     BL_PROFILE("amr-wind::Sampling::fill_buffer");
-    m_scontainer->populate_buffer(m_sample_buf);
+    if (m_out_fmt == "netcdf") {
+#ifdef AMR_WIND_USE_NETCDF
+        m_scontainer->populate_buffer(m_sample_buf);
+#else
+        amrex::Abort(
+            "NetCDF support was not enabled during build time. Please "
+            "recompile or "
+            "use native format");
+#endif
+    }
 }
 
 void Sampling::process_output()
@@ -442,7 +484,7 @@ void Sampling::write_netcdf()
         std::string vname = m_var_names[iv];
         start[1] = 0;
         count[1] = 0;
-        auto offset = iv * num_output_particles();
+        auto offset = iv * num_netcdf_output_particles();
         for (const auto& obj : m_samplers) {
             auto grp = ncf.group(obj->label());
             count[1] = obj->num_output_points();

test/CMakeLists.txt

@@ -263,6 +263,7 @@ if(AMR_WIND_ENABLE_NETCDF)
   add_test_re(abl_meso_input_dpa)
   add_test_re(abl_meso_input_ipa)
   add_test_re(abl_kosovic_neutral)
+  add_test_re(abl_sampling_netcdf)
 endif()
 
 if(AMR_WIND_ENABLE_MASA)

test/test_files/abl_sampling/abl_sampling.inp

@@ -78,6 +78,7 @@ incflo.verbose          =   0          # incflo_level
 #.......................................#
 incflo.post_processing       = sampling sampling2
 sampling.output_frequency    = 1
+sampling.output_format       = native
 sampling.fields              = velocity
 sampling.derived_fields      = mag_vorticity "components(velocity,0,1)" "grad(temperature)" "grad(density)"
 sampling.labels              = volume1
@@ -87,6 +88,7 @@ sampling.volume1.hi          = 500.0 500.0 500.0
 sampling.volume1.num_points  = 30 10 10
 
 sampling2.output_frequency    = 1
+sampling2.output_format       = native
 sampling2.fields              = velocity
 sampling2.derived_fields      = q_criterion_nondim
 sampling2.labels              = volume1

test/test_files/abl_sampling_netcdf/abl_sampling_netcdf.inp

@@ -0,0 +1,98 @@
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            SIMULATION STOP            #
+#.......................................#
+time.stop_time               =   22000.0     # Max (simulated) time to evolve
+time.max_step                =   10          # Max number of time steps
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#         TIME STEP COMPUTATION         #
+#.......................................#
+time.fixed_dt         =   0.5        # Use this constant dt if > 0
+time.cfl              =   0.95         # CFL factor
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            INPUT AND OUTPUT           #
+#.......................................#
+time.plot_interval            =  10       # Steps between plot files
+time.checkpoint_interval      =  5       # Steps between checkpoint files
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               PHYSICS                 #
+#.......................................#
+incflo.gravity          =   0.  0. -9.81  # Gravitational force (3D)
+incflo.density             = 1.0          # Reference density
+
+io.derived_outputs = mag_vorticity "components(velocity,0,1)" "grad(temperature)"
+
+incflo.use_godunov = 1
+incflo.diffusion_type = 2
+transport.viscosity = 1.0e-5
+transport.laminar_prandtl = 0.7
+transport.turbulent_prandtl = 0.3333
+turbulence.model = Smagorinsky
+Smagorinsky_coeffs.Cs = 0.135
+
+
+incflo.physics = ABL
+ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing ABLForcing
+BoussinesqBuoyancy.reference_temperature = 300.0
+ABL.reference_temperature = 300.0
+CoriolisForcing.latitude = 41.3
+ABLForcing.abl_forcing_height = 90
+
+incflo.velocity = 6.128355544951824  5.142300877492314 0.0
+
+ABL.temperature_heights = 650.0 750.0 1000.0
+ABL.temperature_values = 300.0 308.0 308.75
+
+ABL.kappa = .41
+ABL.surface_roughness_z0 = 0.15
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#        ADAPTIVE MESH REFINEMENT       #
+#.......................................#
+amr.n_cell              = 48 48 48    # Grid cells at coarsest AMRlevel
+amr.max_level           = 0           # Max AMR level in hierarchy
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              GEOMETRY                 #
+#.......................................#
+geometry.prob_lo        =   0.       0.     0.  # Lo corner coordinates
+geometry.prob_hi        =   1000.  1000.  1000.  # Hi corner coordinates
+geometry.is_periodic    =   1   1   0   # Periodicity x y z (0/1)
+
+# Boundary conditions
+zlo.type =   "wall_model"
+
+zhi.type =   "slip_wall"
+zhi.temperature_type = "fixed_gradient"
+zhi.temperature = 0.003 # tracer is used to specify potential temperature gradient
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              VERBOSITY                #
+#.......................................#
+incflo.verbose          =   0          # incflo_level
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               SAMPLING                #
+#.......................................#
+incflo.post_processing       = sampling sampling2
+sampling.output_frequency    = 1
+sampling.output_format       = netcdf
+sampling.fields              = velocity
+sampling.derived_fields      = mag_vorticity "components(velocity,0,1)" "grad(temperature)" "grad(density)"
+sampling.labels              = volume1
+sampling.volume1.type        = VolumeSampler
+sampling.volume1.lo          = 200.0 200.0 200.0
+sampling.volume1.hi          = 500.0 500.0 500.0
+sampling.volume1.num_points  = 30 10 10
+
+sampling2.output_frequency    = 1
+sampling.output_format        = netcdf
+sampling2.fields              = velocity
+sampling2.derived_fields      = q_criterion_nondim
+sampling2.labels              = volume1
+sampling2.volume1.type        = VolumeSampler
+sampling2.volume1.lo          = 20.0 20.0 20.0
+sampling2.volume1.hi          = 50.0 50.0 50.0
+sampling2.volume1.num_points  = 5 5 5

tools/fcompare_particles.py

@@ -0,0 +1,51 @@
+import numpy as np
+import argparse
+from pathlib import Path
+import pandas as pd
+import amrex_particle
+from amrex_particle import AmrexParticleFile
+
+
+def main():
+    """Compare particle files"""
+    parser = argparse.ArgumentParser(description="A comparison tool for particles")
+    parser.add_argument("f0", help="A particle directory", type=str)
+    parser.add_argument("f1", help="A particle directory", type=str)
+    parser.add_argument(
+        "-a",
+        "--abs_tol",
+        help="Absolute tolerance (default is 0)",
+        type=float,
+        default=0.0,
+    )
+    parser.add_argument(
+        "-r",
+        "--rel_tol",
+        help="Relative tolerance (default is 0)",
+        type=float,
+        default=0.0,
+    )
+    args = parser.parse_args()
+
+    assert Path(args.f0).is_dir()
+    assert Path(args.f1).is_dir()
+
+    cols_to_drop = ["uid", "set_id", "probe_id"]
+    p0df = AmrexParticleFile(Path(args.f0) / "particles")().drop(cols_to_drop, axis=1)
+    p1df = AmrexParticleFile(Path(args.f1) / "particles")().drop(cols_to_drop, axis=1)
+
+    adiff = np.sqrt(np.square(p0df - p1df).sum(axis=0))
+    rdiff = np.sqrt(np.square(p0df - p1df).sum(axis=0)) / np.sqrt(
+        np.square(p0df).sum(axis=0)
+    )
+    adiff = adiff.to_frame(name="absolute_error")
+    rdiff = rdiff.to_frame(name="relative_error")
+    diff = pd.concat([adiff, rdiff], axis=1).fillna(value={"relative_error": 0.0})
+
+    print(diff)
+    assert (diff["absolute_error"] <= args.abs_tol).all()
+    assert (diff["relative_error"] <= args.rel_tol).all()
+
+
+if __name__ == "__main__":
+    main()