Squared singular values to calculate POD variance (#306)

* Squared singular values

* Added command line option for squared singular values

* Added command line options for nonlinear elasticity

* Fixed typos

* astyle

* comma spacing

* added missing spaces before comma

Author: axla-io

examples/prom/mixed_nonlinear_diffusion.cpp

@@ -391,7 +391,7 @@ void BroadcastUndistributedRomVector(CAROM::Vector* v)
 
 // TODO: move this to the library?
 void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots,
-                std::string name)
+                std::string name, bool squareSV)
 {
     MFEM_VERIFY(nparam > 0, "Must specify a positive number of parameter sets");
 
@@ -411,7 +411,7 @@ void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots,
     generator.endSamples(); // save the merged basis file
 
     int cutoff = 0;
-    generator.finalSummary(1e-4, cutoff, "mergedSV_" + name);
+    generator.finalSummary(1e-4, cutoff, "mergedSV_" + name, 0, squareSV);
 }
 
 std::shared_ptr<const CAROM::Matrix>
@@ -471,6 +471,7 @@ int main(int argc, char *argv[])
     bool use_eqp = false;
     bool writeSampleMesh = false;
     int num_samples_req = -1;
+    bool squareSV = true;
 
     bool pointwiseSnapshots = false;
     int pwx = 0;
@@ -536,6 +537,8 @@ int main(int argc, char *argv[])
                    "Enable or disable the online phase.");
     args.AddOption(&merge, "-merge", "--merge", "-no-merge", "--no-merge",
                    "Enable or disable the merge phase.");
+    args.AddOption(&squareSV, "-sqsv", "--square-sv", "-no-sqsv", "--no-square-sv",
+                   "Use singular values squared in energy fraction.");
     args.AddOption(&samplingType, "-hrtype", "--hrsamplingtype",
                    "Sampling type for hyperreduction.");
     args.AddOption(&num_samples_req, "-nsr", "--nsr",
@@ -664,13 +667,13 @@ int main(int argc, char *argv[])
         totalTimer.Clear();
         totalTimer.Start();
 
-        MergeBasis(R_space.GetTrueVSize(), nsets, max_num_snapshots, "R");
-        MergeBasis(R_space.GetTrueVSize(), nsets, max_num_snapshots, "FR");
-        MergeBasis(W_space.GetTrueVSize(), nsets, max_num_snapshots, "W");
+        MergeBasis(R_space.GetTrueVSize(), nsets, max_num_snapshots, "R", squareSV);
+        MergeBasis(R_space.GetTrueVSize(), nsets, max_num_snapshots, "FR", squareSV);
+        MergeBasis(W_space.GetTrueVSize(), nsets, max_num_snapshots, "W", squareSV);
 
         if (hyperreduce_source)
         {
-            MergeBasis(W_space.GetTrueVSize(), nsets, max_num_snapshots, "S");
+            MergeBasis(W_space.GetTrueVSize(), nsets, max_num_snapshots, "S", squareSV);
         }
 
         totalTimer.Stop();

examples/prom/nonlinear_elasticity_global_rom.cpp

@@ -295,7 +295,7 @@ void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
 
 // TODO: move this to the library?
 void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots,
-                std::string name)
+                std::string name, bool squareSV)
 {
     MFEM_VERIFY(nparam > 0, "Must specify a positive number of parameter sets");
 
@@ -315,7 +315,7 @@ void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots,
     generator.endSamples(); // save the merged basis file
 
     int cutoff = 0;
-    generator.finalSummary(1e-7, cutoff, "mergedSV_" + name + ".txt");
+    generator.finalSummary(1e-7, cutoff, "mergedSV_" + name + ".txt", 0, squareSV);
 }
 
 const CAROM::Matrix *GetSnapshotMatrix(const int dimFOM, const int nparam,
@@ -403,6 +403,7 @@ int main(int argc, char *argv[])
     bool x_base_only = false;
     int num_samples_req = -1;
     const char *samplingType = "gnat";
+    bool squareSV = true;
 
     int nsets = 0;
     int id_param = 0;
@@ -465,6 +466,8 @@ int main(int argc, char *argv[])
                    "Enable or disable the online phase.");
     args.AddOption(&merge, "-merge", "--merge", "-no-merge", "--no-merge",
                    "Enable or disable the merge phase.");
+    args.AddOption(&squareSV, "-sqsv", "--square-sv", "-no-sqsv", "--no-square-sv",
+                   "Use singular values squared in energy fraction.");
     args.AddOption(&samplingType, "-hrtype", "--hrsamplingtype",
                    "Sampling type for hyperreduction.");
     args.AddOption(&num_samples_req, "-nsr", "--nsr",
@@ -630,11 +633,11 @@ int main(int argc, char *argv[])
         // Merge bases
         if (x_base_only == false)
         {
-            MergeBasis(true_size, nsets, max_num_snapshots, "V");
+            MergeBasis(true_size, nsets, max_num_snapshots, "V", squareSV);
         }
 
-        MergeBasis(true_size, nsets, max_num_snapshots, "X");
-        MergeBasis(true_size, nsets, max_num_snapshots, "H");
+        MergeBasis(true_size, nsets, max_num_snapshots, "X", squareSV);
+        MergeBasis(true_size, nsets, max_num_snapshots, "H", squareSV);
 
         totalTimer.Stop();
         if (myid == 0)

lib/linalg/BasisGenerator.cpp

@@ -299,7 +299,7 @@ BasisGenerator::finalSummary(
     const double energyFractionThreshold,
     int & cutoff,
     const std::string & cutoffOutputPath,
-    const int first_sv)
+    const int first_sv, const bool squareSV)
 {
     const int rom_dim = getSpatialBasis()->numColumns();
     std::shared_ptr<const Vector> sing_vals = getSingularValues();
@@ -308,7 +308,8 @@ BasisGenerator::finalSummary(
 
     double sum = 0.0;
     for (int sv = first_sv; sv < sing_vals->dim(); ++sv) {
-        sum += (*sing_vals)(sv);
+        const double s = (*sing_vals)(sv);
+        sum += squareSV ? s * s : s;
     }
 
     int p = std::floor(-std::log10(energyFractionThreshold));
@@ -332,7 +333,8 @@ BasisGenerator::finalSummary(
     }
 
     for (int sv = first_sv; sv < sing_vals->dim(); ++sv) {
-        partialSum += (*sing_vals)(sv);
+        const double s = (*sing_vals)(sv);
+        partialSum += squareSV ? s * s : s;
         for (int i = count; i < p; ++i)
         {
             if (partialSum / sum > 1.0 - std::pow(10, -1 - i))

lib/linalg/BasisGenerator.h

@@ -260,7 +260,7 @@ class BasisGenerator
         const double energyFractionThreshold,
         int & cutoff,
         const std::string & cutoffOutputPath = "",
-        const int first_sv = 0);
+        const int first_sv = 0, const bool squareSV = true);
 
 protected:
     /**