Merge remote-tracking branch 'origin/master' into eigenvalue_example

.github/workflows/run_tests/action.yml

@@ -23,6 +23,10 @@ runs:
           mpirun -n 3 --oversubscribe tests/test_StaticSVD
           ./tests/test_IncrementalSVDBrand
           mpirun -n 3 --oversubscribe tests/test_IncrementalSVDBrand
+          ./tests/test_HDFDatabase
+          mpirun -n 3 --oversubscribe tests/test_HDFDatabase
+          ./tests/test_NNLS
+          mpirun -n 3 --oversubscribe tests/test_NNLS
       shell: bash
     - name: Basis dataset update test
       run: |

CMakeLists.txt

@@ -267,6 +267,7 @@ if(GTEST_FOUND)
   set(unit_test_stems
     Vector
     Matrix
+    HDFDatabase
     DEIM
     DMD
     GNAT
@@ -277,6 +278,7 @@ if(GTEST_FOUND)
     IncrementalSVD
     IncrementalSVDBrand
     GreedyCustomSampler
+    NNLS
     basis_conversion)
   foreach(stem IN LISTS unit_test_stems)
     add_executable(test_${stem} unit_tests/test_${stem}.cpp)

docker/Dockerfile

@@ -10,7 +10,7 @@ WORKDIR /$ENVDIR
 # install packages
 RUN sudo apt-get install -yq git
 RUN sudo apt-get install --no-install-recommends -yq make gcc gfortran libssl-dev cmake
-RUN sudo apt-get install -yq libopenblas-dev libmpich-dev libblas-dev liblapack-dev libscalapack-mpi-dev libhdf5-serial-dev
+RUN sudo apt-get install -yq libopenblas-dev libmpich-dev libblas-dev liblapack-dev libscalapack-mpi-dev libhdf5-mpi-dev hdf5-tools
 RUN sudo apt-get install -yq vim
 RUN sudo apt-get install -yq git-lfs
 RUN sudo apt-get install -yq valgrind

examples/dmd/dg_advection.cpp

@@ -710,9 +710,7 @@ int main(int argc, char *argv[])
     Vector true_solution_u(U->Size());
     true_solution_u = U->GetData();
 
-    dmd_prediction_timer.Start();
-
-    // 14. Predict the state at t_final using DMD.
+    // 14. Predict using DMD.
     if (myid == 0)
     {
         std::cout << "Predicting solution using DMD" << std::endl;
@@ -755,9 +753,9 @@ int main(int argc, char *argv[])
         }
     }
 
-    dmd_prediction_timer.Stop();
-
+    dmd_prediction_timer.Start();
     result_u = dmd_U.predict(t_final);
+    dmd_prediction_timer.Stop();
 
     // 15. Calculate the relative error between the DMD final solution and the true solution.
     Vector dmd_solution_u(result_u->getData(), result_u->dim());

examples/dmd/dg_euler.cpp

@@ -574,9 +574,7 @@ int main(int argc, char *argv[])
     Vector true_solution_e(u_block.GetBlock(3).Size());
     true_solution_e = u_block.GetBlock(3).GetData();
 
-    dmd_prediction_timer.Start();
-
-    // 14. Predict the state at t_final using DMD.
+    // 14. Predict using DMD.
     if (mpi.WorldRank() == 0)
     {
         std::cout << "Predicting density, momentum, and energy using DMD" << std::endl;
@@ -639,12 +637,12 @@ int main(int argc, char *argv[])
         }
     }
 
-    dmd_prediction_timer.Stop();
-
+    dmd_prediction_timer.Start();
     result_dens = dmd_dens->predict(t_final);
     result_x_mom = dmd_x_mom->predict(t_final);
     result_y_mom = dmd_y_mom->predict(t_final);
     result_e = dmd_e->predict(t_final);
+    dmd_prediction_timer.Stop();
 
     // 15. Calculate the relative error between the DMD final solution and the true solution.
     Vector dmd_solution_dens(result_dens->getData(), result_dens->dim());

examples/dmd/heat_conduction.cpp

@@ -498,9 +498,7 @@ int main(int argc, char *argv[])
     Vector true_solution_u(u.Size());
     true_solution_u = u.GetData();
 
-    dmd_prediction_timer.Start();
-
-    // 14. Predict the state at t_final using DMD.
+    // 14. Predict using DMD.
     if (myid == 0)
     {
         std::cout << "Predicting temperature using DMD" << std::endl;
@@ -540,9 +538,9 @@ int main(int argc, char *argv[])
         }
     }
 
-    dmd_prediction_timer.Stop();
-
+    dmd_prediction_timer.Start();
     result_u = dmd_u->predict(t_final);
+    dmd_prediction_timer.Stop();
 
     // 15. Calculate the relative error between the DMD final solution and the true solution.
     Vector dmd_solution_u(result_u->getData(), result_u->dim());

examples/dmd/heat_conduction_dmdc.cpp

@@ -534,9 +534,7 @@ int main(int argc, char *argv[])
     Vector true_solution_u(u.Size());
     true_solution_u = u.GetData();
 
-    dmd_prediction_timer.Start();
-
-    // 14. Predict the state at t_final using DMDc.
+    // 14. Predict using DMDc.
     if (myid == 0)
     {
         std::cout << "Predicting temperature using DMDc" << std::endl;
@@ -576,9 +574,9 @@ int main(int argc, char *argv[])
         }
     }
 
-    dmd_prediction_timer.Stop();
-
+    dmd_prediction_timer.Start();
     result_u = dmd_u->predict(t_final);
+    dmd_prediction_timer.Stop();
 
     // 15. Calculate the relative error between the DMDc final solution and the true solution.
     Vector dmd_solution_u(result_u->getData(), result_u->dim());

examples/dmd/nonlinear_elasticity.cpp

@@ -515,6 +515,12 @@ int main(int argc, char *argv[])
 
         if (last_step || (ti % windowNumSamples) == 0)
         {
+            // Calculate DMD modes
+            if (myid == 0 && rdim != -1 && ef != -1)
+            {
+                std::cout << "Both rdim and ef are set. ef will be ignored." << std::endl;
+            }
+
             if (rdim != -1)
             {
                 if (myid == 0)
@@ -608,10 +614,10 @@ int main(int argc, char *argv[])
         w_gf.Save(ee_ofs);
     }
 
-    // 13. Calculate the DMD modes.
-    if (myid == 0 && rdim != -1 && ef != -1)
+    // 13. Predict using DMD.
+    if (myid == 0)
     {
-        std::cout << "Both rdim and ef are set. ef will be ignored." << std::endl;
+        std::cout << "Predicting position and velocity using DMD" << std::endl;
     }
 
     Vector true_solution_x(x_gf.GetTrueVector().Size());
@@ -620,14 +626,6 @@ int main(int argc, char *argv[])
     Vector true_solution_v(v_gf.GetTrueVector().Size());
     true_solution_v = v_gf.GetTrueVector();
 
-    dmd_prediction_timer.Start();
-
-    // 14. Predict the state at t_final using DMD.
-    if (myid == 0)
-    {
-        std::cout << "Predicting position and velocity using DMD" << std::endl;
-    }
-
     curr_window = 0;
     CAROM::Vector* result_x = dmd_x[curr_window]->predict(ts[0]);
     CAROM::Vector* result_v = dmd_v[curr_window]->predict(ts[0]);
@@ -687,12 +685,12 @@ int main(int argc, char *argv[])
         }
     }
 
-    dmd_prediction_timer.Stop();
-
+    dmd_prediction_timer.Start();
     result_x = dmd_x[curr_window]->predict(t_final);
     result_v = dmd_v[curr_window]->predict(t_final);
+    dmd_prediction_timer.Stop();
 
-    // 15. Calculate the relative error between the DMD final solution and the true solution.
+    // 14. Calculate the relative error between the DMD final solution and the true solution.
     Vector dmd_solution_x(result_x->getData(), result_x->dim());
     Vector diff_x(true_solution_x.Size());
     subtract(dmd_solution_x, true_solution_x, diff_x);
@@ -722,7 +720,7 @@ int main(int argc, char *argv[])
                dmd_prediction_timer.RealTime());
     }
 
-    // 16. Free the used memory.
+    // 15. Free the used memory.
     delete ode_solver;
     delete pmesh;
     delete result_x;

examples/dmd/wave_equation.cpp

@@ -525,8 +525,7 @@ int main(int argc, char *argv[])
         dudt_gf.Save(osol);
     }
 
-    // 10. Predict the state at t_final using DMD.
-    dmd_prediction_timer.Start();
+    // 10. Predict using DMD.
     cout << "Predicting temperature using DMD" << endl;
     CAROM::Vector* result_u = nullptr;
     VisItDataCollection dmd_visit_dc("DMD_Wave_Equation", mesh);
@@ -564,8 +563,9 @@ int main(int argc, char *argv[])
             }
         }
     }
-    dmd_prediction_timer.Stop();
+    dmd_prediction_timer.Start();
     result_u = dmd_u[curr_window]->predict(t_final);
+    dmd_prediction_timer.Stop();
 
     // 11. Calculate the relative error between the DMD final solution and the true solution.
     Vector dmd_solution_u(result_u->getData(), result_u->dim());

examples/misc/combine_samples.cpp

@@ -53,15 +53,19 @@ int main(int argc, char* argv[])
     std::vector<std::string> sample_names;
     int snaps = 0;
     int dim   = 0;
+    int col_min = 1;
+    int col_max = 1e9;
     bool subtract_mean = false;
     bool subtract_offset = false;
+    std::string generator_filename = "total";
     std::string kind = "snapshot";
     std::string offset_file;
 
     bool offset_arg = false;
 
     if (argc >= 2) {
-        for (int i = 1; i < argc; i++) {
+        int i = 1;
+        while (i < argc) {
             if (!strcmp(argv[i], "basis") || !strcmp(argv[i], "-b")) {
                 if (rank==0) std::cout << "Argument " << i << " identified as basis or -b" <<
                                            std::endl;
@@ -80,9 +84,25 @@ int main(int argc, char* argv[])
                 if (rank==0) std::cout << "Will subtract mean" << std::endl;
                 subtract_mean = true;
             }
+            else if (!strcmp(argv[i], "file") || !strcmp(argv[i], "-f")) {
+                i += 1;
+                generator_filename = argv[i];
+                if (rank==0) std::cout << "Output prefix = " << generator_filename << std::endl;
+            }
+            else if (!strcmp(argv[i], "col_min") || !strcmp(argv[i], "-cmin")) {
+                i += 1;
+                col_min = std::stoi(argv[i]);
+                if (rank==0) std::cout << "First column to read = " << col_min << std::endl;
+            }
+            else if (!strcmp(argv[i], "col_max") || !strcmp(argv[i], "-cmax")) {
+                i += 1;
+                col_max = std::stoi(argv[i]);
+                if (rank==0) std::cout << "Last column to read = " << col_max << std::endl;
+            }
             else {
                 sample_names.push_back(argv[i]);
             }
+            i += 1;
         }
     }
     else {
@@ -112,27 +132,29 @@ int main(int argc, char* argv[])
     CAROM_VERIFY((snaps > 0) && (dim > 0));
 
     /*-- Load data from input files --*/
-    std::string generator_filename = "total";
     std::unique_ptr<CAROM::BasisGenerator> static_basis_generator;
     static_basis_generator.reset(new CAROM::BasisGenerator(
                                      CAROM::Options(dim, snaps).setMaxBasisDimension(snaps), false,
                                      generator_filename));
 
     if (rank==0) std::cout << "Loading data from " << kind << std::endl;
     for(const auto& sample_name: sample_names) {
-        static_basis_generator->loadSamples(sample_name, kind);
+        static_basis_generator->loadSampleRange(sample_name, kind, col_min-1,
+                                                col_max-1);
     }
 
     if (rank==0) std::cout << "Saving data uploaded as a snapshot matrix" <<
                                std::endl;
     static_basis_generator->writeSnapshot();
 
+    int rdim;
     if (!subtract_mean && !subtract_offset) {
         /*-- Compute SVD and save file --*/
         if (rank==0) std::cout << "Computing SVD" << std::endl;
         int rom_dim = static_basis_generator->getSpatialBasis()->numColumns();
         if (rank==0) std::cout << "U ROM Dimension: " << rom_dim << std::endl;
         static_basis_generator->endSamples();
+        if (rank==0) static_basis_generator->finalSummary(1e-8, rdim);
     }
     else {
         /*-- load data from hdf5 file to find the mean and subtract it --*/
@@ -188,6 +210,7 @@ int main(int argc, char* argv[])
         int rom_dim = static_basis_generator2->getSpatialBasis()->numColumns();
         if (rank==0) std::cout << "U ROM Dimension: " << rom_dim << std::endl;
         static_basis_generator2->endSamples();
+        if (rank==0) static_basis_generator2->finalSummary(1e-8, rdim);
 
         static_basis_generator2 = nullptr;
     }

examples/prom/dg_advection_global_rom.cpp

@@ -10,14 +10,21 @@
 
 //                       libROM MFEM Example: DG Advection (adapted from ex9p.cpp)
 //
+// =================================================================================
+//
 // Compile with: make dg_advection_global_rom
 //
-// For ROM (reproductive case):
+// Sample runs and results for reproductive case:
 //    dg_advection_global_rom -offline
 //    dg_advection_global_rom -merge -ns 1
 //    dg_advection_global_rom -online
 //
-// For ROM (global rom):
+// Output:
+//  Relative l2 error of ROM solution 0.0005641839849343575
+//
+// =================================================================================
+//
+// Sample runs and results for global ROM:
 // Offline phase: dg_advection_global_rom -offline -ff 1.0 -id 0
 //                dg_advection_global_rom -offline -ff 1.1 -id 1
 //                dg_advection_global_rom -offline -ff 1.2 -id 2
@@ -28,19 +35,10 @@
 //
 // Online phase:  dg_advection_global_rom -online -ff 1.15
 //
-// Sample runs:
-//    mpirun -np 4 dg_advection_global_rom -p 0 -dt 0.005
-//    mpirun -np 4 dg_advection_global_rom -p 0 -dt 0.01
-//    mpirun -np 4 dg_advection_global_rom -p 1 -dt 0.005 -tf 9
-//    mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.002 -tf 9
-//    mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.02 -s 13 -tf 9
-//    mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.004 -tf 9
-//    mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.005 -tf 9
-//    mpirun -np 4 dg_advection_global_rom -p 3 -rp 2 -dt 0.0025 -tf 9 -vs 20
-//    mpirun -np 4 dg_advection_global_rom -p 0 -o 2 -rp 1 -dt 0.01 -tf 8
-//    mpirun -np 4 dg_advection_global_rom -p 0 -rs 2 -dt 0.005 -tf 2
-//    mpirun -np 4 dg_advection_global_rom -p 0 -rs 1 -o 2 -tf 2
-//    mpirun -np 3 dg_advection_global_rom -p 1 -rs 1 -rp 0 -dt 0.005 -tf 0.5
+// Output:
+//   Relative l2 error of ROM solution 0.0004333183604809453
+//
+// =================================================================================
 //
 // Description:  This example code solves the time-dependent advection equation
 //               du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
@@ -128,10 +126,7 @@ class AIR_prec : public Solver
         delete AIR_solver;
         AIR_solver = new HypreBoomerAMG(A_s);
         AIR_solver->SetAdvectiveOptions(1, "", "FA");
-        AIR_solver->SetPrintLevel(
-            0);    // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
-        //    this mesh further in parallel to increase the resolution. Once the
-        //    parallel mesh is defined, the serial mesh can be deleted.
+        AIR_solver->SetPrintLevel(0);
         AIR_solver->SetMaxLevels(50);
     }
 
@@ -464,9 +459,6 @@ int main(int argc, char *argv[])
         {
             cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
         }
-        delete mesh;    // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
-        //    this mesh further in parallel to increase the resolution. Once the
-        //    parallel mesh is defined, the serial mesh can be deleted.
         return 3;
     }