From ff004a498d2b4b384d942937f3e43ce6a9a71613 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 4 Apr 2024 08:51:40 -0700
Subject: [PATCH 01/50] Add boundary condition switch

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index a2087216c..19842870d 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -86,6 +86,7 @@ int main(int argc, char *argv[])
 
     // 2. Parse command-line options.
     const char *mesh_file = "";
+    bool dirichlet = true;
     int ser_ref_levels = 2;
     int par_ref_levels = 1;
     int order = 1;
@@ -117,6 +118,8 @@ int main(int argc, char *argv[])
                    "Mesh file to use.");
     args.AddOption(&problem, "-p", "--problem",
                    "Problem setup to use. See options in Conductivity and Potential functions.");
+    args.AddOption(&dirichlet, "-dir", "--dirichlet", "-neu", "--neumann",
+                   "BC switch.");
     args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
                    "Number of times to refine the mesh uniformly in serial.");
     args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
@@ -349,7 +352,7 @@ int main(int argc, char *argv[])
     if (pmesh->bdr_attributes.Size())
     {
         ess_bdr.SetSize(pmesh->bdr_attributes.Max());
-        ess_bdr = 1;
+        ess_bdr = (dirichlet) ? 1 : 0;
     }
 
     ParBilinearForm *a = new ParBilinearForm(fespace);

From af1f38d950e3ff446ae3f297efb20984f267a22d Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 4 Apr 2024 08:53:56 -0700
Subject: [PATCH 02/50] Output eigenvalue absolute error

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 19842870d..f4c527db1 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -644,12 +644,14 @@ int main(int argc, char *argv[])
         for (int i = 0; i < eigenvalues.Size() && i < nev; i++)
         {
             diff_ev[i] = ev_fom[i] - eigenvalues[i];
-            double ev_diff_norm = sqrt(diff_ev[i] * diff_ev[i]);
-            double ev_fom_norm = sqrt(ev_fom[i] * ev_fom[i]);
             if (myid == 0)
             {
+                std::cout << "FOM solution for eigenvalue " << i << " = " <<
+                          ev_fom[i] << std::endl;
+                std::cout << "Absolute error of ROM solution for eigenvalue " << i << " = " <<
+                          abs(diff_ev[i]) << std::endl;
                 std::cout << "Relative error of ROM solution for eigenvalue " << i << " = " <<
-                          ev_diff_norm / ev_fom_norm << std::endl;
+                          abs(diff_ev[i]) / abs(ev_fom[i]) << std::endl;
             }
         }
 

From 861bcdb3eb9fa6df97db26a4795957bdbf768c9d Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 4 Apr 2024 08:56:09 -0700
Subject: [PATCH 03/50] Make eigenfunction sign in FOM and ROM consistent

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index f4c527db1..29b73f3cc 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -495,6 +495,7 @@ int main(int argc, char *argv[])
     }
 
     // 15. The online phase
+    Vector sign_ev(nev);
     if (online) {
         // 16. read the reduced basis
         assembleTimer.Start();
@@ -677,7 +678,9 @@ int main(int argc, char *argv[])
             mode_fom_ifs.close();
 
             const double fomNorm = sqrt(InnerProduct(mode_fom, mode_fom));
-            mode_fom -= mode_rom;
+            // TODO: use a constant mode_ref for all parameters and (FOM & ROM) eigenfunctions
+            sign_ev[i] = (InnerProduct(mode_fom, mode_rom) >= 0) ? 1 : -1;
+            mode_fom.Add(-sign_ev[i], mode_rom);
             const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom));
             if (myid == 0) std::cout << "Relative l2 error of ROM eigenvector " << i <<
                                          " = " << diffNorm /
@@ -704,6 +707,7 @@ int main(int argc, char *argv[])
                 Vector ev;
                 evect.GetRow(i, ev);
                 x = ev;
+                x *= sign_ev[i];
             }
             visit_evs.push_back(new ParGridFunction(x));
             visit_dc.RegisterField("x" + std::to_string(i), visit_evs.back());
@@ -755,6 +759,7 @@ int main(int argc, char *argv[])
                     Vector ev;
                     evect.GetRow(i, ev);
                     x = ev;
+                    x *= sign_ev[i];
                 }
 
                 sout << "parallel " << num_procs << " " << myid << "\n"

From e60ca5df5685f8b15400352f94c6acadfab7bc7a Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Tue, 23 Apr 2024 12:47:08 -0700
Subject: [PATCH 04/50] Clean up problems

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 157 ++++--------------
 1 file changed, 30 insertions(+), 127 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 29b73f3cc..4b5eb1d7c 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -73,6 +73,8 @@ using namespace mfem;
 double Conductivity(const Vector &x);
 double Potential(const Vector &x);
 int problem = 1;
+int ser_ref_levels = 2;
+int par_ref_levels = 1;
 double kappa;
 Vector bb_min, bb_max;
 
@@ -87,8 +89,6 @@ int main(int argc, char *argv[])
     // 2. Parse command-line options.
     const char *mesh_file = "";
     bool dirichlet = true;
-    int ser_ref_levels = 2;
-    int par_ref_levels = 1;
     int order = 1;
     int nev = 5;
     int seed = 75;
@@ -197,8 +197,6 @@ int main(int argc, char *argv[])
         args.PrintOptions(cout);
     }
 
-    kappa = alpha * M_PI;
-
     if (fom)
     {
         MFEM_VERIFY(fom && !offline && !online
@@ -217,29 +215,14 @@ int main(int argc, char *argv[])
     Mesh *mesh;
     if (mesh_file == "")
     {
-        // square domain with side length 20 at origin (0,0)
-        double x_max = 20.0;
-        double y_max = 20.0;
-        mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false,
-                                              x_max, y_max));
-        mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
-
-        // shift mesh around origin (from [0, bb_max] -> [-bb_max/2, bb_max/2])
-        auto *mesh_transform = +[](const Vector &v_orig, Vector &v_new) -> void {
-            v_new = v_orig;
-            // shift mesh vertices by (bb_max[i] / 2)
-            v_new(0) -= 0.5*bb_max[0];
-            v_new(1) -= 0.5*bb_max[1];
-        };
-
-        // performs the transform - bb_min/bb_max are updated again below
-        mesh->Transform(mesh_transform);
+        mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL));
     }
     else
     {
         mesh = new Mesh(mesh_file, 1, 1);
     }
     int dim = mesh->Dimension();
+    mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
 
     // 4. Refine the mesh in serial to increase the resolution. In this example
     //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
@@ -248,7 +231,6 @@ int main(int argc, char *argv[])
     {
         mesh->UniformRefinement();
     }
-    mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
 
     // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
     //    this mesh further in parallel to increase the resolution. Once the
@@ -338,10 +320,11 @@ int main(int argc, char *argv[])
     assembleTimer.Start();
     ConstantCoefficient one(1.0);
 
+    kappa = alpha * M_PI;
     FunctionCoefficient kappa_0(Conductivity);
     FunctionCoefficient v_0(Potential);
 
-    // Project initial conductivity and potential to visualize initial condition
+    // Project conductivity and potential to visualize initial condition
     ParGridFunction c_gf(fespace);
     c_gf.ProjectCoefficient(kappa_0);
 
@@ -518,23 +501,23 @@ int main(int argc, char *argv[])
                                   numColumnRB);
 
         // 17. form ROM operator
-        CAROM::Matrix invReducedA;
-        ComputeCtAB(*A, *spatialbasis, *spatialbasis, invReducedA);
-        DenseMatrix *A_mat = new DenseMatrix(invReducedA.numRows(),
-                                             invReducedA.numColumns());
-        A_mat->Set(1, invReducedA.getData());
-
-        CAROM::Matrix invReducedM;
-        ComputeCtAB(*M, *spatialbasis, *spatialbasis, invReducedM);
-        DenseMatrix *M_mat = new DenseMatrix(invReducedM.numRows(),
-                                             invReducedM.numColumns());
-        M_mat->Set(1, invReducedM.getData());
+        CAROM::Matrix ReducedA;
+        ComputeCtAB(*A, *spatialbasis, *spatialbasis, ReducedA);
+        DenseMatrix *A_mat = new DenseMatrix(ReducedA.numRows(),
+                                             ReducedA.numColumns());
+        A_mat->Set(1, ReducedA.getData());
+
+        CAROM::Matrix ReducedM;
+        ComputeCtAB(*M, *spatialbasis, *spatialbasis, ReducedM);
+        DenseMatrix *M_mat = new DenseMatrix(ReducedM.numRows(),
+                                             ReducedM.numColumns());
+        M_mat->Set(1, ReducedM.getData());
 
         assembleTimer.Stop();
 
         // 18. solve ROM
         solveTimer.Start();
-        // (Q^T A Q) c = \lamba (Q^T M Q) c
+        // (Q^T A Q) c = \lambda (Q^T M Q) c
         A_mat->Eigenvalues(*M_mat, ev, evect);
         solveTimer.Stop();
 
@@ -694,8 +677,8 @@ int main(int argc, char *argv[])
     VisItDataCollection visit_dc("EllipticEigenproblem", pmesh);
     if (visit)
     {
-        visit_dc.RegisterField("InitialConductivity", &c_gf);
-        visit_dc.RegisterField("InitialPotential", &p_gf);
+        visit_dc.RegisterField("Conductivity", &c_gf);
+        visit_dc.RegisterField("Potential", &p_gf);
         std::vector<ParGridFunction*> visit_evs;
         for (int i = 0; i < nev && i < eigenvalues.Size(); i++)
         {
@@ -837,11 +820,6 @@ double Conductivity(const Vector &x)
         return 1.0 + cos(kappa * x(1)) * sin(kappa * x(0));
     case 3:
     case 4:
-    case 5:
-    case 6:
-    case 7:
-    case 8:
-    case 9:
         return 1.0;
     }
     return 0.0;
@@ -849,101 +827,26 @@ double Conductivity(const Vector &x)
 
 double Potential(const Vector &x)
 {
-    int dim = x.Size();
-
-    Vector X(dim);
-    Vector center(dim); // center of gaussian for problem 4 and 5
-    center = kappa / M_PI;
-
-    Vector neg_center(center);
-    neg_center.Neg();
-
-    // amplitude of gaussians for problems 4-6
-    const double D = 100.0;
-    // width of gaussians for problems 4-6
-    const double c = 0.05;
-
-    const double min_d = 2.0 * sqrt(2.0);
-
-    X = x;
-    for (int i = 0; i < dim; i++)
-    {
-        // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5)
-        center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], center(i));
-        neg_center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i));
-    }
-
+    double D = 100.0;
+    Vector center(x.Size());
     switch (problem)
     {
     case 1:
     case 2:
         return 0.0;
     case 3:
-        return 1.0;
-    case 4:
-        return D * std::exp(-X.DistanceSquaredTo(center) / c);
-    case 5:
-        return -D * std::exp(-X.DistanceSquaredTo(center) / c);
-    case 6:
-        return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp(
-                         -X.DistanceSquaredTo(neg_center) / c));
-    case 7:
-        center = 0.25;
-        center(0) += 0.05 * cos(2.0*kappa);
-        center(1) += 0.05 * sin(2.0*kappa);
-
-        neg_center = center;
-        neg_center.Neg();
-
-        for (int i = 0; i < dim; i++)
+        double c = std::pow(2, -2*(ser_ref_levels + par_ref_levels));
+        for (int i = 0; i < x.Size(); i++)
         {
-            // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5)
-            center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], center(i));
-            neg_center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i));
+            center(i) = 0.5 * (bb_min[i] + bb_max[i]);
         }
-        return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp(
-                         -X.DistanceSquaredTo(neg_center) / c));
-    case 8:
-        // Similar to case 6, but t is restricted to inner 20% of the domain
-        //  in this case, the radius of the gaussian is (0.05*min_d)^2 where
-        //  min_d is the lower bound of the atom distance over time
-
-        // verify t is within inner 20% of domain (centered around mesh origin)
-        CAROM::verify_within_portion(bb_min, bb_max, center, 0.2);
-
-        return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d,
-                              2)) + std::exp(
-                         -X.DistanceSquaredTo(neg_center) / std::pow(c * min_d, 2)));
-    case 9:
-        // Similar to case 7, but t is restricted to inner 20% of the domain and t is defined as:
-        //  t = (1.5 + 0.5*cos(2*k), 1.5 + 0.5*sin(2*k)) where k = alpha * PI, alpha is the input parameter given with the -a option.
-        //  The radius of the gaussian follows case 8: (0.05*min_d)^2
-
-        // Sets the center to vary around +/- 1.5 (absolute location on the mesh)
-        center = 1.5;
-        center(0) += 0.5 * cos(2.0*kappa);
-        center(1) += 0.5 * sin(2.0*kappa);
-
-        // map the absolute location back to a fraction of the mesh domain
-        center(0) = CAROM::map_from_ref_mesh(bb_min[0], bb_max[0], center(0));
-        center(1) = CAROM::map_from_ref_mesh(bb_min[1], bb_max[1], center(1));
-
-        neg_center = center;
-        neg_center.Neg();
-
-        for (int i = 0; i < dim; i++)
+        return D * std::exp(-x.DistanceSquaredTo(center) / c);
+    case 4:
+        for (int i = 0; i < x.Size(); i++)
         {
-            // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5)
-            center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], center(i));
-            neg_center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i));
+            center(i) = 0.5 * (bb_min[i] + bb_max[i]);
         }
-
-        // verify t is within inner 20% of domain (centered around mesh origin)
-        CAROM::verify_within_portion(bb_min, bb_max, center, 0.2);
-
-        return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d,
-                              2)) + std::exp(
-                         -X.DistanceSquaredTo(neg_center) / std::pow(c * min_d, 2)));
+        return D / x.DistanceSquaredTo(center);
     }
     return 0.0;
 }

From c52323708d1b3fad713aefef22b60a062ba2ae89 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Tue, 23 Apr 2024 12:58:30 -0700
Subject: [PATCH 05/50] Change Gaussian radius

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 4b5eb1d7c..7145300fa 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -828,6 +828,7 @@ double Conductivity(const Vector &x)
 double Potential(const Vector &x)
 {
     double D = 100.0;
+    double c = (bb_max[0] - bb_min[0]) * std::pow(2, -2*(ser_ref_levels + par_ref_levels - 1));
     Vector center(x.Size());
     switch (problem)
     {
@@ -835,7 +836,6 @@ double Potential(const Vector &x)
     case 2:
         return 0.0;
     case 3:
-        double c = std::pow(2, -2*(ser_ref_levels + par_ref_levels));
         for (int i = 0; i < x.Size(); i++)
         {
             center(i) = 0.5 * (bb_min[i] + bb_max[i]);

From 333103d5bbd928bbac5b035b6e085eda0731b517 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 24 Apr 2024 07:56:57 -0700
Subject: [PATCH 06/50] Change calculation of Gaussian width

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 7145300fa..2a5cd9577 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -73,9 +73,8 @@ using namespace mfem;
 double Conductivity(const Vector &x);
 double Potential(const Vector &x);
 int problem = 1;
-int ser_ref_levels = 2;
-int par_ref_levels = 1;
 double kappa;
+double gaussian_width;
 Vector bb_min, bb_max;
 
 int main(int argc, char *argv[])
@@ -88,6 +87,8 @@ int main(int argc, char *argv[])
 
     // 2. Parse command-line options.
     const char *mesh_file = "";
+    int ser_ref_levels = 2;
+    int par_ref_levels = 1;
     bool dirichlet = true;
     int order = 1;
     int nev = 5;
@@ -223,6 +224,7 @@ int main(int argc, char *argv[])
     }
     int dim = mesh->Dimension();
     mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
+    gaussian_width = mesh->GetElementVolume(0) * pow(2.0, 1.0 - ser_ref_levels - par_ref_levels);
 
     // 4. Refine the mesh in serial to increase the resolution. In this example
     //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
@@ -828,7 +830,6 @@ double Conductivity(const Vector &x)
 double Potential(const Vector &x)
 {
     double D = 100.0;
-    double c = (bb_max[0] - bb_min[0]) * std::pow(2, -2*(ser_ref_levels + par_ref_levels - 1));
     Vector center(x.Size());
     switch (problem)
     {
@@ -840,7 +841,7 @@ double Potential(const Vector &x)
         {
             center(i) = 0.5 * (bb_min[i] + bb_max[i]);
         }
-        return D * std::exp(-x.DistanceSquaredTo(center) / c);
+        return D * std::exp(-x.DistanceSquaredTo(center) / pow(gaussian_width, 2.0));
     case 4:
         for (int i = 0; i < x.Size(); i++)
         {

From 0e1f90383759d2352fd7e212f60f6d7830317ff9 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 22 May 2024 14:15:57 -0700
Subject: [PATCH 07/50] Slight changes in example command lines

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 2a5cd9577..7eb9c0a07 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -32,7 +32,7 @@
 // Merge phase:   elliptic_eigenproblem_global_rom -merge -p 2 -ns 4
 //
 // FOM run (for error calculation):
-//                elliptic_eigenproblem_global_rom -fom -p 2 -a 0.5
+//                elliptic_eigenproblem_global_rom -fom -p 2 -a 0.50
 // Example output:
 //   Number of unknowns: 289
 //    Eigenvalue 0: 0.04533314
@@ -43,7 +43,7 @@
 //   Elapsed time for assembling FOM: 1.471708e-03 second
 //   Elapsed time for solving FOM: 3.416310e-01 second
 //
-// Online phase:  elliptic_eigenproblem_global_rom -online -p 2 -a 0.5
+// Online phase:  elliptic_eigenproblem_global_rom -online -p 2 -a 0.50
 // Example output:
 //   Eigenvalue 0: = 0.048430949
 //   Eigenvalue 1: = 0.12021157

From 8239b14b4490d27499a571735107dc7fc30393bf Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 22 May 2024 14:18:31 -0700
Subject: [PATCH 08/50] Change filenames

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 7eb9c0a07..e731c6fd6 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -271,8 +271,8 @@ int main(int argc, char *argv[])
     int max_num_snapshots = 100;
     bool update_right_SV = false;
     bool isIncremental = false;
-    const std::string basisName = "basis";
-    const std::string basisFileName = basisName + std::to_string(id);
+    const std::string baseName = "elliptic_eigenproblem_";
+    const std::string basis_filename = baseName + "basis";
     CAROM::Options* options;
     CAROM::BasisGenerator *generator;
     StopWatch solveTimer, assembleTimer, mergeTimer;
@@ -282,7 +282,8 @@ int main(int argc, char *argv[])
     {
         options = new CAROM::Options(fespace->GetTrueVSize(), nev, nev,
                                      update_right_SV);
-        generator = new CAROM::BasisGenerator(*options, isIncremental, basisFileName);
+        std::string snapshot_basename = baseName + "par" + std::to_string(id);
+        generator = new CAROM::BasisGenerator(*options, isIncremental, snapshot_basename);
     }
 
     // 9. The merge phase
@@ -291,10 +292,10 @@ int main(int argc, char *argv[])
         mergeTimer.Start();
         options = new CAROM::Options(fespace->GetTrueVSize(), max_num_snapshots, nev,
                                      update_right_SV);
-        generator = new CAROM::BasisGenerator(*options, isIncremental, basisName);
+        generator = new CAROM::BasisGenerator(*options, isIncremental, basis_filename);
         for (int paramID=0; paramID<nsets; ++paramID)
         {
-            std::string snapshot_filename = basisName + std::to_string(
+            std::string snapshot_filename = baseName + "par" + std::to_string(
                                                 paramID) + "_snapshot";
             generator->loadSamples(snapshot_filename, "snapshot");
         }
@@ -484,7 +485,7 @@ int main(int argc, char *argv[])
     if (online) {
         // 16. read the reduced basis
         assembleTimer.Start();
-        CAROM::BasisReader reader(basisName);
+        CAROM::BasisReader reader(basis_filename);
 
         Vector ev;
         const CAROM::Matrix *spatialbasis;

From 157a2e99d66f58adfecdf36020515ea77c59ad6a Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 22 May 2024 14:20:44 -0700
Subject: [PATCH 09/50] Remove TODO

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 1 -
 1 file changed, 1 deletion(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index e731c6fd6..956c35389 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -664,7 +664,6 @@ int main(int argc, char *argv[])
             mode_fom_ifs.close();
 
             const double fomNorm = sqrt(InnerProduct(mode_fom, mode_fom));
-            // TODO: use a constant mode_ref for all parameters and (FOM & ROM) eigenfunctions
             sign_ev[i] = (InnerProduct(mode_fom, mode_rom) >= 0) ? 1 : -1;
             mode_fom.Add(-sign_ev[i], mode_rom);
             const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom));

From 6c7a0f8ba4c4bd9b6726fe2da502ac5b2ddeb952 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 22 May 2024 14:23:05 -0700
Subject: [PATCH 10/50] Astyle

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 956c35389..f94dd44d5 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -224,7 +224,8 @@ int main(int argc, char *argv[])
     }
     int dim = mesh->Dimension();
     mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
-    gaussian_width = mesh->GetElementVolume(0) * pow(2.0, 1.0 - ser_ref_levels - par_ref_levels);
+    gaussian_width = mesh->GetElementVolume(0) * pow(2.0,
+                     1.0 - ser_ref_levels - par_ref_levels);
 
     // 4. Refine the mesh in serial to increase the resolution. In this example
     //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
@@ -283,7 +284,8 @@ int main(int argc, char *argv[])
         options = new CAROM::Options(fespace->GetTrueVSize(), nev, nev,
                                      update_right_SV);
         std::string snapshot_basename = baseName + "par" + std::to_string(id);
-        generator = new CAROM::BasisGenerator(*options, isIncremental, snapshot_basename);
+        generator = new CAROM::BasisGenerator(*options, isIncremental,
+                                              snapshot_basename);
     }
 
     // 9. The merge phase

From 8ab3481a7cf81618ac149caab426aa706acfb9ac Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 29 May 2024 22:37:32 -0700
Subject: [PATCH 11/50] Change filenames. Change calculation of Gaussian width

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 24 ++++++++++++-------
 1 file changed, 15 insertions(+), 9 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index f94dd44d5..9c9ff7bc9 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -224,8 +224,9 @@ int main(int argc, char *argv[])
     }
     int dim = mesh->Dimension();
     mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
-    gaussian_width = mesh->GetElementVolume(0) * pow(2.0,
-                     1.0 - ser_ref_levels - par_ref_levels);
+    double h_min, h_max, k_min, k_max;
+    mesh->GetCharacteristics(h_min, h_max, k_min, k_max);
+    gaussian_width = h_max * pow(2.0, 1.0 - ser_ref_levels - par_ref_levels);
 
     // 4. Refine the mesh in serial to increase the resolution. In this example
     //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
@@ -637,6 +638,8 @@ int main(int argc, char *argv[])
             {
                 std::cout << "FOM solution for eigenvalue " << i << " = " <<
                           ev_fom[i] << std::endl;
+                std::cout << "ROM solution for eigenvalue " << i << " = " <<
+                          eigenvalues[i] << std::endl;
                 std::cout << "Absolute error of ROM solution for eigenvalue " << i << " = " <<
                           abs(diff_ev[i]) << std::endl;
                 std::cout << "Relative error of ROM solution for eigenvalue " << i << " = " <<
@@ -678,11 +681,14 @@ int main(int argc, char *argv[])
         }
     }
 
-    VisItDataCollection visit_dc("EllipticEigenproblem", pmesh);
+    VisItDataCollection *visit_dc = NULL;
     if (visit)
     {
-        visit_dc.RegisterField("Conductivity", &c_gf);
-        visit_dc.RegisterField("Potential", &p_gf);
+        if (offline) visit_dc = new VisItDataCollection("elliptic_eigenproblem_offline_par" + std::to_string(id), pmesh);
+        else if (fom) visit_dc = new VisItDataCollection("elliptic_eigenproblem_fom", pmesh);
+        else if (online) visit_dc = new VisItDataCollection("elliptic_eigenproblem_rom", pmesh);
+        visit_dc->RegisterField("Conductivity", &c_gf);
+        visit_dc->RegisterField("Potential", &p_gf);
         std::vector<ParGridFunction*> visit_evs;
         for (int i = 0; i < nev && i < eigenvalues.Size(); i++)
         {
@@ -697,11 +703,11 @@ int main(int argc, char *argv[])
                 x *= sign_ev[i];
             }
             visit_evs.push_back(new ParGridFunction(x));
-            visit_dc.RegisterField("x" + std::to_string(i), visit_evs.back());
+            visit_dc->RegisterField("Eigenmode_" + std::to_string(i), visit_evs.back());
         }
-        visit_dc.SetCycle(0);
-        visit_dc.SetTime(0.0);
-        visit_dc.Save();
+        visit_dc->SetCycle(0);
+        visit_dc->SetTime(0.0);
+        visit_dc->Save();
         for (size_t i = 0; i < visit_evs.size(); i++)
         {
             delete visit_evs[i];

From e2add1b0ef9706f5223992ee5554437dabdad920 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 30 May 2024 08:53:59 -0700
Subject: [PATCH 12/50] Slight changes to filenames

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 9c9ff7bc9..d9dacc2f5 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -684,9 +684,9 @@ int main(int argc, char *argv[])
     VisItDataCollection *visit_dc = NULL;
     if (visit)
     {
-        if (offline) visit_dc = new VisItDataCollection("elliptic_eigenproblem_offline_par" + std::to_string(id), pmesh);
-        else if (fom) visit_dc = new VisItDataCollection("elliptic_eigenproblem_fom", pmesh);
-        else if (online) visit_dc = new VisItDataCollection("elliptic_eigenproblem_rom", pmesh);
+        if (offline) visit_dc = new VisItDataCollection(baseName + "offline_par" + std::to_string(id), pmesh);
+        else if (fom) visit_dc = new VisItDataCollection(baseName + "fom", pmesh);
+        else if (online) visit_dc = new VisItDataCollection(baseName + "rom", pmesh);
         visit_dc->RegisterField("Conductivity", &c_gf);
         visit_dc->RegisterField("Potential", &p_gf);
         std::vector<ParGridFunction*> visit_evs;

From 8e17778f008d80987e67ce9081bd5b94cc70e62e Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 30 May 2024 09:35:28 -0700
Subject: [PATCH 13/50] Consistent sign over parameters and eigenvector
 ordering

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 62 ++++++++++++++-----
 1 file changed, 48 insertions(+), 14 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index d9dacc2f5..ed0bc6b34 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -484,7 +484,6 @@ int main(int argc, char *argv[])
     }
 
     // 15. The online phase
-    Vector sign_ev(nev);
     if (online) {
         // 16. read the reduced basis
         assembleTimer.Start();
@@ -567,8 +566,9 @@ int main(int argc, char *argv[])
 
     // 19. Save the refined mesh and the modes in parallel. This output can be
     //     viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
+    Vector sign_ev(nev);
     {
-        ostringstream mesh_name, mode_name;
+        ostringstream mesh_name, mode_name, mode_ref_name;
         mesh_name << "elliptic_eigenproblem-mesh." << setfill('0') << setw(6) << myid;
 
         ofstream mesh_ofs(mesh_name.str().c_str());
@@ -576,26 +576,60 @@ int main(int argc, char *argv[])
         pmesh->Print(mesh_ofs);
 
         std::string mode_prefix = "mode_";
-        if (online)
-        {
-            mode_prefix += "rom_";
-        } else if (fom)
+        if (fom || offline)
         {
             mode_prefix += "fom_";
         }
+        else if (online)
+        {
+            mode_prefix += "rom_";
+        }
 
+        Vector mode_ref(evect.NumCols());
         for (int i=0; i < nev && i < eigenvalues.Size(); i++)
         {
+            if (!(offline && id == 0))
+            {
+                mode_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
+                      << setfill('0') << setw(6) << myid;
+
+                ifstream mode_ref_ifs(mode_name.str().c_str());
+                mode_ref_ifs.precision(16);
+                mode_ref.Load(mode_ref_ifs, evect.NumCols());
+                mode_ref_ifs.close();
+            }
+
+            Vector ev;
             if (fom || offline) {
-                // convert eigenvector from HypreParVector to ParGridFunction
-                x = lobpcg->GetEigenvector(i);
-            } else {
+                ev = lobpcg->GetEigenvector(i);
+            }
+            else
+            {
                 // for online, eigenvectors are stored in evect matrix
-                Vector ev;
                 evect.GetRow(i, ev);
-                x = ev;
             }
+            sign_ev[i] = (InnerProduct(mode_ref, ev) >= 0) ? 1 : -1;
+            ev *= sign_ev[i];
+            // convert eigenvector from HypreParVector to ParGridFunction
+            x = ev;
 
+            if (offline && id == 0)
+            {
+                mode_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
+                          << setfill('0') << setw(6) << myid;
+
+                ofstream mode_ofs(mode_name.str().c_str());
+                mode_ofs.precision(16);
+
+                // TODO: issue using .Load() function if file written with .Save()?
+                //x.Save(mode_ofs);
+                for (int j = 0; j < x.Size(); j++)
+                {
+                    mode_ofs << x[j] << "\n";
+                }
+                mode_ofs.close();
+            }
+ 
             mode_name << mode_prefix << setfill('0') << setw(2) << i << "."
                       << setfill('0') << setw(6) << myid;
 
@@ -669,8 +703,7 @@ int main(int argc, char *argv[])
             mode_fom_ifs.close();
 
             const double fomNorm = sqrt(InnerProduct(mode_fom, mode_fom));
-            sign_ev[i] = (InnerProduct(mode_fom, mode_rom) >= 0) ? 1 : -1;
-            mode_fom.Add(-sign_ev[i], mode_rom);
+            mode_fom.Add(-1.0, mode_rom);
             const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom));
             if (myid == 0) std::cout << "Relative l2 error of ROM eigenvector " << i <<
                                          " = " << diffNorm /
@@ -695,7 +728,8 @@ int main(int argc, char *argv[])
             if (fom || offline) {
                 // convert eigenvector from HypreParVector to ParGridFunction
                 x = lobpcg->GetEigenvector(i);
-            } else {
+            }
+            else {
                 // for online, eigenvectors are stored in evect matrix
                 Vector ev;
                 evect.GetRow(i, ev);

From 69812cfc5fbaa622d0ddc0f3e5115d39d539782b Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 30 May 2024 09:44:57 -0700
Subject: [PATCH 14/50] Astyle

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index ed0bc6b34..7bf8deaf3 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -591,7 +591,7 @@ int main(int argc, char *argv[])
             if (!(offline && id == 0))
             {
                 mode_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
-                      << setfill('0') << setw(6) << myid;
+                          << setfill('0') << setw(6) << myid;
 
                 ifstream mode_ref_ifs(mode_name.str().c_str());
                 mode_ref_ifs.precision(16);
@@ -629,7 +629,7 @@ int main(int argc, char *argv[])
                 }
                 mode_ofs.close();
             }
- 
+
             mode_name << mode_prefix << setfill('0') << setw(2) << i << "."
                       << setfill('0') << setw(6) << myid;
 
@@ -717,7 +717,8 @@ int main(int argc, char *argv[])
     VisItDataCollection *visit_dc = NULL;
     if (visit)
     {
-        if (offline) visit_dc = new VisItDataCollection(baseName + "offline_par" + std::to_string(id), pmesh);
+        if (offline) visit_dc = new VisItDataCollection(baseName + "offline_par" +
+                    std::to_string(id), pmesh);
         else if (fom) visit_dc = new VisItDataCollection(baseName + "fom", pmesh);
         else if (online) visit_dc = new VisItDataCollection(baseName + "rom", pmesh);
         visit_dc->RegisterField("Conductivity", &c_gf);

From 958bd256c3e6c98d17b0f997a5c73564018ac32b Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 30 May 2024 10:17:13 -0700
Subject: [PATCH 15/50] Fix visualization signs

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 25 ++++++++++---------
 1 file changed, 13 insertions(+), 12 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 7bf8deaf3..7ee2bd19c 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -726,17 +726,18 @@ int main(int argc, char *argv[])
         std::vector<ParGridFunction*> visit_evs;
         for (int i = 0; i < nev && i < eigenvalues.Size(); i++)
         {
+            Vector ev;
             if (fom || offline) {
-                // convert eigenvector from HypreParVector to ParGridFunction
-                x = lobpcg->GetEigenvector(i);
+                ev = lobpcg->GetEigenvector(i);
             }
             else {
                 // for online, eigenvectors are stored in evect matrix
-                Vector ev;
                 evect.GetRow(i, ev);
-                x = ev;
-                x *= sign_ev[i];
             }
+            // convert eigenvector from HypreParVector to ParGridFunction
+            ev *= sign_ev[i];
+            x = ev;
+            // convert eigenvector from HypreParVector to ParGridFunction
             visit_evs.push_back(new ParGridFunction(x));
             visit_dc->RegisterField("Eigenmode_" + std::to_string(i), visit_evs.back());
         }
@@ -779,16 +780,16 @@ int main(int argc, char *argv[])
                          << ", Lambda = " << eigenvalues[i] << endl;
                 }
 
+                Vector ev;
                 if (fom || offline) {
-                    // convert eigenvector from HypreParVector to ParGridFunction
-                    x = lobpcg->GetEigenvector(i);
-                } else {
-                    // for online, eigenvectors are stored in evect matrix
-                    Vector ev;
+                    ev = lobpcg->GetEigenvector(i);
+                }
+                else {
                     evect.GetRow(i, ev);
-                    x = ev;
-                    x *= sign_ev[i];
                 }
+                // convert eigenvector from HypreParVector to ParGridFunction
+                ev *= sign_ev[i];
+                x = ev;
 
                 sout << "parallel " << num_procs << " " << myid << "\n"
                      << "solution\n" << *pmesh << x << flush

From ec0d4c1cbceda21bd08605c23a58bd977b968902 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 13 Jun 2024 08:40:13 -0700
Subject: [PATCH 16/50] Fix stream names

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 20 +++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 7ee2bd19c..e1582762f 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -590,10 +590,10 @@ int main(int argc, char *argv[])
         {
             if (!(offline && id == 0))
             {
-                mode_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
-                          << setfill('0') << setw(6) << myid;
+                mode_ref_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
+                              << setfill('0') << setw(6) << myid;
 
-                ifstream mode_ref_ifs(mode_name.str().c_str());
+                ifstream mode_ref_ifs(mode_ref_name.str().c_str());
                 mode_ref_ifs.precision(16);
                 mode_ref.Load(mode_ref_ifs, evect.NumCols());
                 mode_ref_ifs.close();
@@ -615,19 +615,19 @@ int main(int argc, char *argv[])
 
             if (offline && id == 0)
             {
-                mode_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
-                          << setfill('0') << setw(6) << myid;
+                mode_ref_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
+                              << setfill('0') << setw(6) << myid;
 
-                ofstream mode_ofs(mode_name.str().c_str());
-                mode_ofs.precision(16);
+                ofstream mode_ref_ofs(mode_ref_name.str().c_str());
+                mode_ref_ofs.precision(16);
 
                 // TODO: issue using .Load() function if file written with .Save()?
-                //x.Save(mode_ofs);
+                //x.Save(mode_ref_ofs);
                 for (int j = 0; j < x.Size(); j++)
                 {
-                    mode_ofs << x[j] << "\n";
+                    mode_ref_ofs << x[j] << "\n";
                 }
-                mode_ofs.close();
+                mode_ref_ofs.close();
             }
 
             mode_name << mode_prefix << setfill('0') << setw(2) << i << "."

From 9e7805253453a95ea9caecae351d9f50edc1fadd Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 13 Jun 2024 08:47:56 -0700
Subject: [PATCH 17/50] Reset stream name

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index e1582762f..2a508df65 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -590,6 +590,7 @@ int main(int argc, char *argv[])
         {
             if (!(offline && id == 0))
             {
+                mode_ref_name.str("");
                 mode_ref_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
                               << setfill('0') << setw(6) << myid;
 
@@ -615,6 +616,7 @@ int main(int argc, char *argv[])
 
             if (offline && id == 0)
             {
+                mode_ref_name.str("");
                 mode_ref_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
                               << setfill('0') << setw(6) << myid;
 

From 1e864d6ef3470dfaeabdc9b97f5227b037386383 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 13 Jun 2024 09:23:19 -0700
Subject: [PATCH 18/50] Reorganize reference mode sign logic

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 43 +++++++++----------
 1 file changed, 21 insertions(+), 22 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 2a508df65..bea059285 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -585,23 +585,11 @@ int main(int argc, char *argv[])
             mode_prefix += "rom_";
         }
 
-        Vector mode_ref(evect.NumCols());
         for (int i=0; i < nev && i < eigenvalues.Size(); i++)
         {
-            if (!(offline && id == 0))
-            {
-                mode_ref_name.str("");
-                mode_ref_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
-                              << setfill('0') << setw(6) << myid;
-
-                ifstream mode_ref_ifs(mode_ref_name.str().c_str());
-                mode_ref_ifs.precision(16);
-                mode_ref.Load(mode_ref_ifs, evect.NumCols());
-                mode_ref_ifs.close();
-            }
-
             Vector ev;
-            if (fom || offline) {
+            if (fom || offline)
+            {
                 ev = lobpcg->GetEigenvector(i);
             }
             else
@@ -609,17 +597,13 @@ int main(int argc, char *argv[])
                 // for online, eigenvectors are stored in evect matrix
                 evect.GetRow(i, ev);
             }
-            sign_ev[i] = (InnerProduct(mode_ref, ev) >= 0) ? 1 : -1;
-            ev *= sign_ev[i];
-            // convert eigenvector from HypreParVector to ParGridFunction
-            x = ev;
+            Vector mode_ref(ev.Size());
 
+            mode_ref_name.str("");
+            mode_ref_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
+                              << setfill('0') << setw(6) << myid;
             if (offline && id == 0)
             {
-                mode_ref_name.str("");
-                mode_ref_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
-                              << setfill('0') << setw(6) << myid;
-
                 ofstream mode_ref_ofs(mode_ref_name.str().c_str());
                 mode_ref_ofs.precision(16);
 
@@ -631,7 +615,22 @@ int main(int argc, char *argv[])
                 }
                 mode_ref_ofs.close();
             }
+            else
+            {
+                ifstream mode_ref_ifs(mode_ref_name.str().c_str());
+                mode_ref_ifs.precision(16);
+                mode_ref.Load(mode_ref_ifs, evect.NumCols());
+                mode_ref_ifs.close();
+            }
+            sign_ev[i] = (InnerProduct(mode_ref, ev) >= 0) ? 1 : -1;
+            std::cout << "norm_ref[i] = " << InnerProduct(mode_ref, mode_ref) << std::endl;
+            std::cout << "norm_ev[i] = " << InnerProduct(ev, ev) << std::endl;
+            std::cout << "ip[i] = " << InnerProduct(mode_ref, ev) << std::endl;
+            ev *= sign_ev[i];
+            // convert eigenvector from HypreParVector to ParGridFunction
+            x = ev;
 
+            mode_name.str("");
             mode_name << mode_prefix << setfill('0') << setw(2) << i << "."
                       << setfill('0') << setw(6) << myid;
 

From 95d9533f636fb56c1e80f7e8796f197cf44ac273 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 13 Jun 2024 09:27:42 -0700
Subject: [PATCH 19/50] Clean up

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index bea059285..b15eeb8f8 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -623,9 +623,6 @@ int main(int argc, char *argv[])
                 mode_ref_ifs.close();
             }
             sign_ev[i] = (InnerProduct(mode_ref, ev) >= 0) ? 1 : -1;
-            std::cout << "norm_ref[i] = " << InnerProduct(mode_ref, mode_ref) << std::endl;
-            std::cout << "norm_ev[i] = " << InnerProduct(ev, ev) << std::endl;
-            std::cout << "ip[i] = " << InnerProduct(mode_ref, ev) << std::endl;
             ev *= sign_ev[i];
             // convert eigenvector from HypreParVector to ParGridFunction
             x = ev;

From 8760daf9cdadfc806285db5e922ebca1a6b65e9a Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 13 Jun 2024 12:43:33 -0700
Subject: [PATCH 20/50] Change variable names

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 43 ++++++++++---------
 1 file changed, 23 insertions(+), 20 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index b15eeb8f8..50849fc25 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -14,9 +14,8 @@
 //
 // Description:  This example code demonstrates the use of MFEM and libROM to
 //               define a simple projection-based reduced order model of the
-//               eigenvalue problem -Delta ((1 + alpha v) u) = lambda u with homogeneous
-//               Dirichlet boundary conditions where alpha is a scalar ROM parameter
-//               controlling the frequency of v.
+//               eigenvalue problem Lu = lambda u with homogeneous
+//               Dirichlet boundary conditions. 
 //
 //               We compute a number of the lowest eigenmodes by discretizing
 //               the Laplacian and Mass operators using a FE space of the
@@ -73,7 +72,7 @@ using namespace mfem;
 double Conductivity(const Vector &x);
 double Potential(const Vector &x);
 int problem = 1;
-double kappa;
+double gaussian_depth = 1.0;
 double gaussian_width;
 Vector bb_min, bb_max;
 
@@ -97,7 +96,6 @@ int main(int argc, char *argv[])
     bool sp_solver = false;
     bool cpardiso_solver = false;
 
-    double alpha = 1.0;
     bool visualization = true;
     bool visit = false;
     int vis_steps = 5;
@@ -131,8 +129,8 @@ int main(int argc, char *argv[])
                    "Number of desired eigenmodes.");
     args.AddOption(&seed, "-s", "--seed",
                    "Random seed used to initialize LOBPCG.");
-    args.AddOption(&alpha, "-a", "--alpha",
-                   "Alpha coefficient.");
+    args.AddOption(&gaussian_depth, "-a", "--gaussian-depth",
+                   "Gaussian depth.");
     args.AddOption(&id, "-id", "--id", "Parametric id");
     args.AddOption(&nsets, "-ns", "--nset", "Number of parametric snapshot sets");
     args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -326,7 +324,6 @@ int main(int argc, char *argv[])
     assembleTimer.Start();
     ConstantCoefficient one(1.0);
 
-    kappa = alpha * M_PI;
     FunctionCoefficient kappa_0(Conductivity);
     FunctionCoefficient v_0(Potential);
 
@@ -597,13 +594,14 @@ int main(int argc, char *argv[])
                 // for online, eigenvectors are stored in evect matrix
                 evect.GetRow(i, ev);
             }
-            Vector mode_ref(ev.Size());
 
-            mode_ref_name.str("");
+            std::cout << "ev.Size = " << ev.Size() << std::endl;
+            Vector mode_ref(ev.Size());
             mode_ref_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
-                              << setfill('0') << setw(6) << myid;
+                          << setfill('0') << setw(6) << myid;
             if (offline && id == 0)
             {
+                mode_ref = ev;
                 ofstream mode_ref_ofs(mode_ref_name.str().c_str());
                 mode_ref_ofs.precision(16);
 
@@ -611,7 +609,7 @@ int main(int argc, char *argv[])
                 //x.Save(mode_ref_ofs);
                 for (int j = 0; j < x.Size(); j++)
                 {
-                    mode_ref_ofs << x[j] << "\n";
+                    mode_ref_ofs << mode_ref[j] << "\n";
                 }
                 mode_ref_ofs.close();
             }
@@ -619,15 +617,18 @@ int main(int argc, char *argv[])
             {
                 ifstream mode_ref_ifs(mode_ref_name.str().c_str());
                 mode_ref_ifs.precision(16);
-                mode_ref.Load(mode_ref_ifs, evect.NumCols());
+                mode_ref.Load(mode_ref_ifs, ev.Size());
                 mode_ref_ifs.close();
             }
+            mode_ref_name.str("");
+            std::cout << i << "-th norm ref = " << InnerProduct(mode_ref, mode_ref) << std::endl;
+            std::cout << i << "-th norm ev = " << InnerProduct(ev, ev) << std::endl;
+            std::cout << i << "-th inner product = " << InnerProduct(mode_ref, ev) << std::endl;
             sign_ev[i] = (InnerProduct(mode_ref, ev) >= 0) ? 1 : -1;
             ev *= sign_ev[i];
             // convert eigenvector from HypreParVector to ParGridFunction
             x = ev;
 
-            mode_name.str("");
             mode_name << mode_prefix << setfill('0') << setw(2) << i << "."
                       << setfill('0') << setw(6) << myid;
 
@@ -854,14 +855,17 @@ int main(int argc, char *argv[])
 
 double Conductivity(const Vector &x)
 {
-    int dim = x.Size();
-
+    Vector center(x.Size());
     switch (problem)
     {
     case 1:
         return 1.0;
     case 2:
-        return 1.0 + cos(kappa * x(1)) * sin(kappa * x(0));
+        for (int i = 0; i < x.Size(); i++)
+        {
+            center(i) = 0.5 * (bb_min[i] + bb_max[i]);
+        }
+        return 1.0 + gaussian_depth * std::exp(-x.DistanceSquaredTo(center) / pow(gaussian_width, 2.0));
     case 3:
     case 4:
         return 1.0;
@@ -871,7 +875,6 @@ double Conductivity(const Vector &x)
 
 double Potential(const Vector &x)
 {
-    double D = 100.0;
     Vector center(x.Size());
     switch (problem)
     {
@@ -883,13 +886,13 @@ double Potential(const Vector &x)
         {
             center(i) = 0.5 * (bb_min[i] + bb_max[i]);
         }
-        return D * std::exp(-x.DistanceSquaredTo(center) / pow(gaussian_width, 2.0));
+        return gaussian_depth * std::exp(-x.DistanceSquaredTo(center) / pow(gaussian_width, 2.0));
     case 4:
         for (int i = 0; i < x.Size(); i++)
         {
             center(i) = 0.5 * (bb_min[i] + bb_max[i]);
         }
-        return D / x.DistanceSquaredTo(center);
+        return gaussian_depth / x.DistanceSquaredTo(center);
     }
     return 0.0;
 }

From 1ebd1f280bcc8525aac06aed789c7b26fccf8f36 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 13 Jun 2024 13:37:12 -0700
Subject: [PATCH 21/50] Warning meesages

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 25 ++++++++++++-------
 1 file changed, 16 insertions(+), 9 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 50849fc25..9a3df093d 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -90,7 +90,7 @@ int main(int argc, char *argv[])
     int par_ref_levels = 1;
     bool dirichlet = true;
     int order = 1;
-    int nev = 5;
+    int nev = 4;
     int seed = 75;
     bool slu_solver = false;
     bool sp_solver = false;
@@ -573,13 +573,16 @@ int main(int argc, char *argv[])
         pmesh->Print(mesh_ofs);
 
         std::string mode_prefix = "mode_";
+        std::string mode_ref_prefix = "mode_";
         if (fom || offline)
         {
             mode_prefix += "fom_";
+            mode_ref_prefix += "ref_";
         }
         else if (online)
         {
             mode_prefix += "rom_";
+            mode_ref_prefix += "fom_";
         }
 
         for (int i=0; i < nev && i < eigenvalues.Size(); i++)
@@ -595,9 +598,8 @@ int main(int argc, char *argv[])
                 evect.GetRow(i, ev);
             }
 
-            std::cout << "ev.Size = " << ev.Size() << std::endl;
             Vector mode_ref(ev.Size());
-            mode_ref_name << "mode_ref_" << setfill('0') << setw(2) << i << "."
+            mode_ref_name << mode_ref_prefix << setfill('0') << setw(2) << i << "."
                           << setfill('0') << setw(6) << myid;
             if (offline && id == 0)
             {
@@ -621,9 +623,14 @@ int main(int argc, char *argv[])
                 mode_ref_ifs.close();
             }
             mode_ref_name.str("");
-            std::cout << i << "-th norm ref = " << InnerProduct(mode_ref, mode_ref) << std::endl;
-            std::cout << i << "-th norm ev = " << InnerProduct(ev, ev) << std::endl;
-            std::cout << i << "-th inner product = " << InnerProduct(mode_ref, ev) << std::endl;
+            if (abs(InnerProduct(mode_ref, ev)) < 0.9 * InnerProduct(mode_ref, mode_ref)) 
+            {
+                std::cout << "Warning: " << i << "-th eigenvector is not directly comparable." << std::endl;
+                std::cout << "TODO: Projection error of eigenvector." << std::endl;
+                std::cout << "Square norm of reference " << i << " = " << InnerProduct(mode_ref, mode_ref) << std::endl;
+                std::cout << "Square norm of eigenvector " << i << " = " << InnerProduct(ev, ev) << std::endl;
+                std::cout << "Inner product = " << InnerProduct(mode_ref, ev) << std::endl;
+            }
             sign_ev[i] = (InnerProduct(mode_ref, ev) >= 0) ? 1 : -1;
             ev *= sign_ev[i];
             // convert eigenvector from HypreParVector to ParGridFunction
@@ -726,17 +733,17 @@ int main(int argc, char *argv[])
         for (int i = 0; i < nev && i < eigenvalues.Size(); i++)
         {
             Vector ev;
-            if (fom || offline) {
+            if (fom || offline)
+            {
                 ev = lobpcg->GetEigenvector(i);
             }
             else {
                 // for online, eigenvectors are stored in evect matrix
                 evect.GetRow(i, ev);
             }
-            // convert eigenvector from HypreParVector to ParGridFunction
             ev *= sign_ev[i];
-            x = ev;
             // convert eigenvector from HypreParVector to ParGridFunction
+            x = ev;
             visit_evs.push_back(new ParGridFunction(x));
             visit_dc->RegisterField("Eigenmode_" + std::to_string(i), visit_evs.back());
         }

From 91fba19b278ac9be87a98869eafd25e3bd675ac6 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 13 Jun 2024 13:38:02 -0700
Subject: [PATCH 22/50] Astyle

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 19 ++++++++++++-------
 1 file changed, 12 insertions(+), 7 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 9a3df093d..46f1e6908 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -15,7 +15,7 @@
 // Description:  This example code demonstrates the use of MFEM and libROM to
 //               define a simple projection-based reduced order model of the
 //               eigenvalue problem Lu = lambda u with homogeneous
-//               Dirichlet boundary conditions. 
+//               Dirichlet boundary conditions.
 //
 //               We compute a number of the lowest eigenmodes by discretizing
 //               the Laplacian and Mass operators using a FE space of the
@@ -623,12 +623,15 @@ int main(int argc, char *argv[])
                 mode_ref_ifs.close();
             }
             mode_ref_name.str("");
-            if (abs(InnerProduct(mode_ref, ev)) < 0.9 * InnerProduct(mode_ref, mode_ref)) 
+            if (abs(InnerProduct(mode_ref, ev)) < 0.9 * InnerProduct(mode_ref, mode_ref))
             {
-                std::cout << "Warning: " << i << "-th eigenvector is not directly comparable." << std::endl;
+                std::cout << "Warning: " << i << "-th eigenvector is not directly comparable."
+                          << std::endl;
                 std::cout << "TODO: Projection error of eigenvector." << std::endl;
-                std::cout << "Square norm of reference " << i << " = " << InnerProduct(mode_ref, mode_ref) << std::endl;
-                std::cout << "Square norm of eigenvector " << i << " = " << InnerProduct(ev, ev) << std::endl;
+                std::cout << "Square norm of reference " << i << " = " << InnerProduct(mode_ref,
+                          mode_ref) << std::endl;
+                std::cout << "Square norm of eigenvector " << i << " = " << InnerProduct(ev,
+                          ev) << std::endl;
                 std::cout << "Inner product = " << InnerProduct(mode_ref, ev) << std::endl;
             }
             sign_ev[i] = (InnerProduct(mode_ref, ev) >= 0) ? 1 : -1;
@@ -872,7 +875,8 @@ double Conductivity(const Vector &x)
         {
             center(i) = 0.5 * (bb_min[i] + bb_max[i]);
         }
-        return 1.0 + gaussian_depth * std::exp(-x.DistanceSquaredTo(center) / pow(gaussian_width, 2.0));
+        return 1.0 + gaussian_depth * std::exp(-x.DistanceSquaredTo(center) / pow(
+                gaussian_width, 2.0));
     case 3:
     case 4:
         return 1.0;
@@ -893,7 +897,8 @@ double Potential(const Vector &x)
         {
             center(i) = 0.5 * (bb_min[i] + bb_max[i]);
         }
-        return gaussian_depth * std::exp(-x.DistanceSquaredTo(center) / pow(gaussian_width, 2.0));
+        return gaussian_depth * std::exp(-x.DistanceSquaredTo(center) / pow(
+                                             gaussian_width, 2.0));
     case 4:
         for (int i = 0; i < x.Size(); i++)
         {

From 11d63b0b512e52bc5546969fd4c13af3bb85fa88 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 13 Jun 2024 15:04:14 -0700
Subject: [PATCH 23/50] Normalize

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 26 ++++++++-----------
 1 file changed, 11 insertions(+), 15 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 46f1e6908..0320a9c7f 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -623,28 +623,25 @@ int main(int argc, char *argv[])
                 mode_ref_ifs.close();
             }
             mode_ref_name.str("");
-            if (abs(InnerProduct(mode_ref, ev)) < 0.9 * InnerProduct(mode_ref, mode_ref))
+            sign_ev[i] = (InnerProduct(mode_ref, ev) >= 0) ? 1 : -1;
+            ev *= sign_ev[i] / sqrt(InnerProduct(ev, ev));
+
+            if (InnerProduct(mode_ref, ev) < 0.9)
             {
-                std::cout << "Warning: " << i << "-th eigenvector is not directly comparable."
+                std::cout << "Warning: eigenvector " << i <<
+                          " in FOM and ROM are not directly comparable."
                           << std::endl;
-                std::cout << "TODO: Projection error of eigenvector." << std::endl;
-                std::cout << "Square norm of reference " << i << " = " << InnerProduct(mode_ref,
-                          mode_ref) << std::endl;
-                std::cout << "Square norm of eigenvector " << i << " = " << InnerProduct(ev,
-                          ev) << std::endl;
                 std::cout << "Inner product = " << InnerProduct(mode_ref, ev) << std::endl;
+                std::cout << "TODO: Projection error of eigenvector." << std::endl;
             }
-            sign_ev[i] = (InnerProduct(mode_ref, ev) >= 0) ? 1 : -1;
-            ev *= sign_ev[i];
+
             // convert eigenvector from HypreParVector to ParGridFunction
             x = ev;
-
             mode_name << mode_prefix << setfill('0') << setw(2) << i << "."
                       << setfill('0') << setw(6) << myid;
 
             ofstream mode_ofs(mode_name.str().c_str());
             mode_ofs.precision(16);
-
             // TODO: issue using .Load() function if file written with .Save()?
             //x.Save(mode_ofs);
             for (int j = 0; j < x.Size(); j++)
@@ -715,8 +712,7 @@ int main(int argc, char *argv[])
             mode_fom.Add(-1.0, mode_rom);
             const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom));
             if (myid == 0) std::cout << "Relative l2 error of ROM eigenvector " << i <<
-                                         " = " << diffNorm /
-                                         fomNorm << std::endl;
+                                         " = " << diffNorm / fomNorm << std::endl;
 
             mode_name.str("");
             mode_name_fom.str("");
@@ -744,7 +740,7 @@ int main(int argc, char *argv[])
                 // for online, eigenvectors are stored in evect matrix
                 evect.GetRow(i, ev);
             }
-            ev *= sign_ev[i];
+            ev *= sign_ev[i] / sqrt(InnerProduct(ev, ev));
             // convert eigenvector from HypreParVector to ParGridFunction
             x = ev;
             visit_evs.push_back(new ParGridFunction(x));
@@ -797,7 +793,7 @@ int main(int argc, char *argv[])
                     evect.GetRow(i, ev);
                 }
                 // convert eigenvector from HypreParVector to ParGridFunction
-                ev *= sign_ev[i];
+                ev *= sign_ev[i] / sqrt(InnerProduct(ev, ev));
                 x = ev;
 
                 sout << "parallel " << num_procs << " " << myid << "\n"

From a819f828b4b6188d64bb52b42e8d42a37a865f80 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 19 Jun 2024 12:58:51 -0700
Subject: [PATCH 24/50] Adjust Gaussian

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 88 +++++++++++++------
 1 file changed, 62 insertions(+), 26 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 0320a9c7f..34af63c92 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -72,9 +72,12 @@ using namespace mfem;
 double Conductivity(const Vector &x);
 double Potential(const Vector &x);
 int problem = 1;
-double gaussian_depth = 1.0;
+double gaussian_depth = -800.0;
+double gaussian_center = 0.0;
 double gaussian_width;
+Vector center;
 Vector bb_min, bb_max;
+double h_min, h_max, k_min, k_max;
 
 int main(int argc, char *argv[])
 {
@@ -131,6 +134,8 @@ int main(int argc, char *argv[])
                    "Random seed used to initialize LOBPCG.");
     args.AddOption(&gaussian_depth, "-a", "--gaussian-depth",
                    "Gaussian depth.");
+    args.AddOption(&gaussian_center, "-c", "--gaussian-center",
+                   "Gaussian center.");
     args.AddOption(&id, "-id", "--id", "Parametric id");
     args.AddOption(&nsets, "-ns", "--nset", "Number of parametric snapshot sets");
     args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -222,9 +227,29 @@ int main(int argc, char *argv[])
     }
     int dim = mesh->Dimension();
     mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
-    double h_min, h_max, k_min, k_max;
     mesh->GetCharacteristics(h_min, h_max, k_min, k_max);
-    gaussian_width = h_max * pow(2.0, 1.0 - ser_ref_levels - par_ref_levels);
+    h_max *= pow(0.5, ser_ref_levels + par_ref_levels);
+
+    gaussian_width = 5.0 * h_max;
+    center.SetSize(dim);
+    center(0) = h_max * gaussian_center + 0.5 * (bb_min[0] + bb_max[0]);
+    for (int i = 1; i < dim; i++)
+    {
+        center(i) = 0.5 * (bb_min[i] + bb_max[i]);
+    }
+    if (myid == 0)
+    {
+        std::cout << "Bounding box = [" << bb_min[0] << ", " << bb_max[0] << "]";
+        for (int i = 1; i < dim; i++)
+            std::cout << " x [" << bb_min[i] << ", " << bb_max[i] << "]";
+        std::cout << std::endl << "Mesh size = " << h_max << std::endl;
+
+        std::cout << "Gaussian depth = " << gaussian_depth
+                  << ", width = " << gaussian_width
+                  << ", center = (" << center(0);
+        for (int i = 1; i < dim; i++) std::cout << ", " << center(i);
+        std::cout << ")" << std::endl;
+    }
 
     // 4. Refine the mesh in serial to increase the resolution. In this example
     //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
@@ -481,7 +506,8 @@ int main(int argc, char *argv[])
     }
 
     // 15. The online phase
-    if (online) {
+    if (online)
+    {
         // 16. read the reduced basis
         assembleTimer.Start();
         CAROM::BasisReader reader(basis_filename);
@@ -693,28 +719,51 @@ int main(int argc, char *argv[])
         Vector mode_fom(evect.NumCols());
         for (int i = 0; i < eigenvalues.Size() && i < nev; i++)
         {
-            mode_name << "mode_rom_" << setfill('0') << setw(2) << i << "."
-                      << setfill('0') << setw(6) << myid;
             mode_name_fom << "mode_fom_" << setfill('0') << setw(2) << i << "."
                           << setfill('0') << setw(6) << myid;
 
-            ifstream mode_rom_ifs(mode_name.str().c_str());
-            mode_rom_ifs.precision(16);
-            mode_rom.Load(mode_rom_ifs, evect.NumCols());
-            mode_rom_ifs.close();
-
             ifstream mode_fom_ifs(mode_name_fom.str().c_str());
             mode_fom_ifs.precision(16);
             mode_fom.Load(mode_fom_ifs, evect.NumCols());
             mode_fom_ifs.close();
 
             const double fomNorm = sqrt(InnerProduct(mode_fom, mode_fom));
-            mode_fom.Add(-1.0, mode_rom);
+
+            for (int j = 0; j < eigenvalues.Size() && j < nev; j++)
+            {
+                if (abs(ev_fom[j] - ev_fom[i]) < 1e-6)
+                {
+                    mode_name << "mode_rom_" << setfill('0') << setw(2) << j << "."
+                              << setfill('0') << setw(6) << myid;
+                    ifstream mode_rom_ifs(mode_name.str().c_str());
+                    mode_rom_ifs.precision(16);
+                    mode_rom.Load(mode_rom_ifs, evect.NumCols());
+                    mode_rom_ifs.close();
+                    std::cout << "Projecting FOM mode " << i 
+                              << " onto ROM mode " << j << std::endl;
+                    std::cout << "ip = " <<  InnerProduct(mode_fom, mode_rom) << std::endl;
+                    mode_fom.Add(-InnerProduct(mode_fom, mode_rom), mode_rom);
+                    mode_name.str("");
+                }
+                else
+                {
+                    mode_name << "mode_rom_" << setfill('0') << setw(2) << j << "."
+                              << setfill('0') << setw(6) << myid;
+                    ifstream mode_rom_ifs(mode_name.str().c_str());
+                    mode_rom_ifs.precision(16);
+                    mode_rom.Load(mode_rom_ifs, evect.NumCols());
+                    mode_rom_ifs.close();
+                    std::cout << "Not projecting FOM mode " << i 
+                              << " onto ROM mode " << j << std::endl;
+                    std::cout << "ip = " <<  InnerProduct(mode_fom, mode_rom) << std::endl;
+                    mode_name.str("");
+                }
+            }
+
             const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom));
             if (myid == 0) std::cout << "Relative l2 error of ROM eigenvector " << i <<
                                          " = " << diffNorm / fomNorm << std::endl;
 
-            mode_name.str("");
             mode_name_fom.str("");
         }
     }
@@ -867,10 +916,6 @@ double Conductivity(const Vector &x)
     case 1:
         return 1.0;
     case 2:
-        for (int i = 0; i < x.Size(); i++)
-        {
-            center(i) = 0.5 * (bb_min[i] + bb_max[i]);
-        }
         return 1.0 + gaussian_depth * std::exp(-x.DistanceSquaredTo(center) / pow(
                 gaussian_width, 2.0));
     case 3:
@@ -882,24 +927,15 @@ double Conductivity(const Vector &x)
 
 double Potential(const Vector &x)
 {
-    Vector center(x.Size());
     switch (problem)
     {
     case 1:
     case 2:
         return 0.0;
     case 3:
-        for (int i = 0; i < x.Size(); i++)
-        {
-            center(i) = 0.5 * (bb_min[i] + bb_max[i]);
-        }
         return gaussian_depth * std::exp(-x.DistanceSquaredTo(center) / pow(
                                              gaussian_width, 2.0));
     case 4:
-        for (int i = 0; i < x.Size(); i++)
-        {
-            center(i) = 0.5 * (bb_min[i] + bb_max[i]);
-        }
         return gaussian_depth / x.DistanceSquaredTo(center);
     }
     return 0.0;

From b4fc1fd094eca82239840d37e763a4921ab1d339 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 19 Jun 2024 15:22:18 -0700
Subject: [PATCH 25/50] Rename variables

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 198 +++++++++---------
 1 file changed, 99 insertions(+), 99 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 34af63c92..ab32723e6 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -72,9 +72,8 @@ using namespace mfem;
 double Conductivity(const Vector &x);
 double Potential(const Vector &x);
 int problem = 1;
-double gaussian_depth = -800.0;
-double gaussian_center = 0.0;
-double gaussian_width;
+double amplitude = -800.0;
+double relative_center = 0.0;
 Vector center;
 Vector bb_min, bb_max;
 double h_min, h_max, k_min, k_max;
@@ -132,10 +131,10 @@ int main(int argc, char *argv[])
                    "Number of desired eigenmodes.");
     args.AddOption(&seed, "-s", "--seed",
                    "Random seed used to initialize LOBPCG.");
-    args.AddOption(&gaussian_depth, "-a", "--gaussian-depth",
-                   "Gaussian depth.");
-    args.AddOption(&gaussian_center, "-c", "--gaussian-center",
-                   "Gaussian center.");
+    args.AddOption(&amplitude, "-a", "--amplitude",
+                   "Amplitude of coefficient fields.");
+    args.AddOption(&relative_center, "-c", "--center",
+                   "Number of grid elements to which center is shifted.");
     args.AddOption(&id, "-id", "--id", "Parametric id");
     args.AddOption(&nsets, "-ns", "--nset", "Number of parametric snapshot sets");
     args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -230,25 +229,10 @@ int main(int argc, char *argv[])
     mesh->GetCharacteristics(h_min, h_max, k_min, k_max);
     h_max *= pow(0.5, ser_ref_levels + par_ref_levels);
 
-    gaussian_width = 5.0 * h_max;
     center.SetSize(dim);
-    center(0) = h_max * gaussian_center + 0.5 * (bb_min[0] + bb_max[0]);
-    for (int i = 1; i < dim; i++)
+    for (int i = 0; i < dim; i++)
     {
-        center(i) = 0.5 * (bb_min[i] + bb_max[i]);
-    }
-    if (myid == 0)
-    {
-        std::cout << "Bounding box = [" << bb_min[0] << ", " << bb_max[0] << "]";
-        for (int i = 1; i < dim; i++)
-            std::cout << " x [" << bb_min[i] << ", " << bb_max[i] << "]";
-        std::cout << std::endl << "Mesh size = " << h_max << std::endl;
-
-        std::cout << "Gaussian depth = " << gaussian_depth
-                  << ", width = " << gaussian_width
-                  << ", center = (" << center(0);
-        for (int i = 1; i < dim; i++) std::cout << ", " << center(i);
-        std::cout << ")" << std::endl;
+        center(i) = h_max * relative_center + 0.5 * (bb_min[i] + bb_max[i]);
     }
 
     // 4. Refine the mesh in serial to increase the resolution. In this example
@@ -407,7 +391,7 @@ int main(int argc, char *argv[])
     ParGridFunction x(fespace);
     HypreLOBPCG *lobpcg;
     Array<double> eigenvalues;
-    DenseMatrix evect;
+    DenseMatrix eigenvectors;
 
     // 11. The offline phase
     if(fom || offline)
@@ -512,7 +496,7 @@ int main(int argc, char *argv[])
         assembleTimer.Start();
         CAROM::BasisReader reader(basis_filename);
 
-        Vector ev;
+        Vector eigenvalues_rom;
         const CAROM::Matrix *spatialbasis;
         if (rdim != -1)
         {
@@ -546,28 +530,30 @@ int main(int argc, char *argv[])
         // 18. solve ROM
         solveTimer.Start();
         // (Q^T A Q) c = \lambda (Q^T M Q) c
-        A_mat->Eigenvalues(*M_mat, ev, evect);
+        A_mat->Eigenvalues(*M_mat, eigenvalues_rom, eigenvectors);
         solveTimer.Stop();
 
         if (myid == 0)
         {
-            eigenvalues = Array<double>(ev.GetData(), ev.Size());
-            for (int i = 0; i < ev.Size() && i < nev; i++)
+            eigenvalues = Array<double>(eigenvalues_rom.GetData(), eigenvalues_rom.Size());
+            for (int i = 0; i < eigenvalues_rom.Size() && i < nev; i++)
             {
                 std::cout << "Eigenvalue " << i << ": = " << eigenvalues[i] << "\n";
             }
         }
 
-        DenseMatrix tmp(evect);
-        evect = DenseMatrix(nev, numRowRB);
-        for (int i = 0; i < ev.Size() && i < nev; i++)
+        DenseMatrix tmp(eigenvectors);
+        eigenvectors = DenseMatrix(nev, numRowRB);
+        for (int i = 0; i < eigenvalues_rom.Size() && i < nev; i++)
         {
-            Vector evector;
-            tmp.GetRow(i, evector);
-            CAROM::Vector evector_carom(evector.GetData(), evector.Size(), false, false);
-            CAROM::Vector *ev_carom = spatialbasis->mult(evector_carom);
-            evect.SetRow(i, ev_carom->getData());
-            delete ev_carom;
+            Vector reduced_eigenvector;
+            tmp.GetRow(i, reduced_eigenvector);
+            CAROM::Vector reduced_eigenvector_carom(reduced_eigenvector.GetData(),
+                                                    reduced_eigenvector.Size(), false, false);
+            CAROM::Vector *eigenvector_carom = spatialbasis->mult(
+                                                   reduced_eigenvector_carom);
+            eigenvectors.SetRow(i, eigenvector_carom->getData());
+            delete eigenvector_carom;
         }
 
         delete spatialbasis;
@@ -576,20 +562,22 @@ int main(int argc, char *argv[])
         delete M_mat;
     }
 
-    ostringstream sol_ev_name, sol_ev_name_fom;
+    ostringstream sol_eigenvalue_name, sol_eigenvalue_name_fom;
     if (fom || offline)
     {
-        sol_ev_name << "sol_eigenvalues_fom." << setfill('0') << setw(6) << myid;
+        sol_eigenvalue_name << "sol_eigenvalues_fom." << setfill('0') << setw(
+                                6) << myid;
     }
     if (online)
     {
-        sol_ev_name << "sol_eigenvalues." << setfill('0') << setw(6) << myid;
-        sol_ev_name_fom << "sol_eigenvalues_fom." << setfill('0') << setw(6) << myid;
+        sol_eigenvalue_name << "sol_eigenvalues." << setfill('0') << setw(6) << myid;
+        sol_eigenvalue_name_fom << "sol_eigenvalues_fom." << setfill('0') << setw(
+                                    6) << myid;
     }
 
     // 19. Save the refined mesh and the modes in parallel. This output can be
     //     viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
-    Vector sign_ev(nev);
+    Vector sign_eigenvectors(nev);
     {
         ostringstream mesh_name, mode_name, mode_ref_name;
         mesh_name << "elliptic_eigenproblem-mesh." << setfill('0') << setw(6) << myid;
@@ -613,23 +601,23 @@ int main(int argc, char *argv[])
 
         for (int i=0; i < nev && i < eigenvalues.Size(); i++)
         {
-            Vector ev;
+            Vector eigenvector_i;
             if (fom || offline)
             {
-                ev = lobpcg->GetEigenvector(i);
+                eigenvector_i = lobpcg->GetEigenvector(i);
             }
             else
             {
-                // for online, eigenvectors are stored in evect matrix
-                evect.GetRow(i, ev);
+                // for online, eigenvectors are stored in eigenvectors matrix
+                eigenvectors.GetRow(i, eigenvector_i);
             }
 
-            Vector mode_ref(ev.Size());
+            Vector mode_ref(eigenvector_i.Size());
             mode_ref_name << mode_ref_prefix << setfill('0') << setw(2) << i << "."
                           << setfill('0') << setw(6) << myid;
             if (offline && id == 0)
             {
-                mode_ref = ev;
+                mode_ref = eigenvector_i;
                 ofstream mode_ref_ofs(mode_ref_name.str().c_str());
                 mode_ref_ofs.precision(16);
 
@@ -645,24 +633,26 @@ int main(int argc, char *argv[])
             {
                 ifstream mode_ref_ifs(mode_ref_name.str().c_str());
                 mode_ref_ifs.precision(16);
-                mode_ref.Load(mode_ref_ifs, ev.Size());
+                mode_ref.Load(mode_ref_ifs, eigenvector_i.Size());
                 mode_ref_ifs.close();
             }
             mode_ref_name.str("");
-            sign_ev[i] = (InnerProduct(mode_ref, ev) >= 0) ? 1 : -1;
-            ev *= sign_ev[i] / sqrt(InnerProduct(ev, ev));
+            sign_eigenvectors[i] = (InnerProduct(mode_ref, eigenvector_i) >= 0) ? 1 : -1;
+            eigenvector_i *= sign_eigenvectors[i] / sqrt(InnerProduct(eigenvector_i,
+                             eigenvector_i));
 
-            if (InnerProduct(mode_ref, ev) < 0.9)
+            if (InnerProduct(mode_ref, eigenvector_i) < 0.9)
             {
                 std::cout << "Warning: eigenvector " << i <<
                           " in FOM and ROM are not directly comparable."
                           << std::endl;
-                std::cout << "Inner product = " << InnerProduct(mode_ref, ev) << std::endl;
-                std::cout << "TODO: Projection error of eigenvector." << std::endl;
+                std::cout << "Inner product = " << InnerProduct(mode_ref,
+                          eigenvector_i) << std::endl;
+                std::cout << "TODO: Visualization of projected eigenvector." << std::endl;
             }
 
             // convert eigenvector from HypreParVector to ParGridFunction
-            x = ev;
+            x = eigenvector_i;
             mode_name << mode_prefix << setfill('0') << setw(2) << i << "."
                       << setfill('0') << setw(6) << myid;
 
@@ -678,45 +668,45 @@ int main(int argc, char *argv[])
             mode_name.str("");
         }
 
-        ofstream sol_ev_ofs(sol_ev_name.str().c_str());
-        sol_ev_ofs.precision(16);
+        ofstream sol_eigenvalue_ofs(sol_eigenvalue_name.str().c_str());
+        sol_eigenvalue_ofs.precision(16);
         for (int i = 0; i < nev && i < eigenvalues.Size(); ++i)
         {
-            sol_ev_ofs << eigenvalues[i] << std::endl;
+            sol_eigenvalue_ofs << eigenvalues[i] << std::endl;
         }
     }
 
     if (online)
     {
         // Initialize FOM solution
-        Vector ev_fom(nev);
+        Vector eigenvalues_fom(nev);
 
         ifstream fom_file;
-        fom_file.open(sol_ev_name_fom.str().c_str());
-        ev_fom.Load(fom_file, ev_fom.Size());
+        fom_file.open(sol_eigenvalue_name_fom.str().c_str());
+        eigenvalues_fom.Load(fom_file, eigenvalues_fom.Size());
         fom_file.close();
 
-        Vector diff_ev(nev);
+        Vector diff_eigenvalues(nev);
         for (int i = 0; i < eigenvalues.Size() && i < nev; i++)
         {
-            diff_ev[i] = ev_fom[i] - eigenvalues[i];
+            diff_eigenvalues[i] = eigenvalues_fom[i] - eigenvalues[i];
             if (myid == 0)
             {
                 std::cout << "FOM solution for eigenvalue " << i << " = " <<
-                          ev_fom[i] << std::endl;
+                          eigenvalues_fom[i] << std::endl;
                 std::cout << "ROM solution for eigenvalue " << i << " = " <<
                           eigenvalues[i] << std::endl;
                 std::cout << "Absolute error of ROM solution for eigenvalue " << i << " = " <<
-                          abs(diff_ev[i]) << std::endl;
+                          abs(diff_eigenvalues[i]) << std::endl;
                 std::cout << "Relative error of ROM solution for eigenvalue " << i << " = " <<
-                          abs(diff_ev[i]) / abs(ev_fom[i]) << std::endl;
+                          abs(diff_eigenvalues[i]) / abs(eigenvalues_fom[i]) << std::endl;
             }
         }
 
         // Calculate errors of eigenvectors
         ostringstream mode_name, mode_name_fom;
-        Vector mode_rom(evect.NumCols());
-        Vector mode_fom(evect.NumCols());
+        Vector mode_rom(eigenvectors.NumCols());
+        Vector mode_fom(eigenvectors.NumCols());
         for (int i = 0; i < eigenvalues.Size() && i < nev; i++)
         {
             mode_name_fom << "mode_fom_" << setfill('0') << setw(2) << i << "."
@@ -724,22 +714,22 @@ int main(int argc, char *argv[])
 
             ifstream mode_fom_ifs(mode_name_fom.str().c_str());
             mode_fom_ifs.precision(16);
-            mode_fom.Load(mode_fom_ifs, evect.NumCols());
+            mode_fom.Load(mode_fom_ifs, eigenvectors.NumCols());
             mode_fom_ifs.close();
 
             const double fomNorm = sqrt(InnerProduct(mode_fom, mode_fom));
 
             for (int j = 0; j < eigenvalues.Size() && j < nev; j++)
             {
-                if (abs(ev_fom[j] - ev_fom[i]) < 1e-6)
+                if (abs(eigenvalues_fom[j] - eigenvalues_fom[i]) < 1e-6)
                 {
                     mode_name << "mode_rom_" << setfill('0') << setw(2) << j << "."
                               << setfill('0') << setw(6) << myid;
                     ifstream mode_rom_ifs(mode_name.str().c_str());
                     mode_rom_ifs.precision(16);
-                    mode_rom.Load(mode_rom_ifs, evect.NumCols());
+                    mode_rom.Load(mode_rom_ifs, eigenvectors.NumCols());
                     mode_rom_ifs.close();
-                    std::cout << "Projecting FOM mode " << i 
+                    std::cout << "Projecting FOM mode " << i
                               << " onto ROM mode " << j << std::endl;
                     std::cout << "ip = " <<  InnerProduct(mode_fom, mode_rom) << std::endl;
                     mode_fom.Add(-InnerProduct(mode_fom, mode_rom), mode_rom);
@@ -751,9 +741,9 @@ int main(int argc, char *argv[])
                               << setfill('0') << setw(6) << myid;
                     ifstream mode_rom_ifs(mode_name.str().c_str());
                     mode_rom_ifs.precision(16);
-                    mode_rom.Load(mode_rom_ifs, evect.NumCols());
+                    mode_rom.Load(mode_rom_ifs, eigenvectors.NumCols());
                     mode_rom_ifs.close();
-                    std::cout << "Not projecting FOM mode " << i 
+                    std::cout << "Not projecting FOM mode " << i
                               << " onto ROM mode " << j << std::endl;
                     std::cout << "ip = " <<  InnerProduct(mode_fom, mode_rom) << std::endl;
                     mode_name.str("");
@@ -777,30 +767,32 @@ int main(int argc, char *argv[])
         else if (online) visit_dc = new VisItDataCollection(baseName + "rom", pmesh);
         visit_dc->RegisterField("Conductivity", &c_gf);
         visit_dc->RegisterField("Potential", &p_gf);
-        std::vector<ParGridFunction*> visit_evs;
+        std::vector<ParGridFunction*> visit_eigenvectors;
         for (int i = 0; i < nev && i < eigenvalues.Size(); i++)
         {
-            Vector ev;
+            Vector eigenvector_i;
             if (fom || offline)
             {
-                ev = lobpcg->GetEigenvector(i);
+                eigenvector_i = lobpcg->GetEigenvector(i);
             }
             else {
-                // for online, eigenvectors are stored in evect matrix
-                evect.GetRow(i, ev);
+                // for online, eigenvectors are stored in eigvenvectors matrix
+                eigenvectors.GetRow(i, eigenvector_i);
             }
-            ev *= sign_ev[i] / sqrt(InnerProduct(ev, ev));
+            eigenvector_i *= sign_eigenvectors[i] / sqrt(InnerProduct(eigenvector_i,
+                             eigenvector_i));
             // convert eigenvector from HypreParVector to ParGridFunction
-            x = ev;
-            visit_evs.push_back(new ParGridFunction(x));
-            visit_dc->RegisterField("Eigenmode_" + std::to_string(i), visit_evs.back());
+            x = eigenvector_i;
+            visit_eigenvectors.push_back(new ParGridFunction(x));
+            visit_dc->RegisterField("Eigenmode_" + std::to_string(i),
+                                    visit_eigenvectors.back());
         }
         visit_dc->SetCycle(0);
         visit_dc->SetTime(0.0);
         visit_dc->Save();
-        for (size_t i = 0; i < visit_evs.size(); i++)
+        for (size_t i = 0; i < visit_eigenvectors.size(); i++)
         {
-            delete visit_evs[i];
+            delete visit_eigenvectors[i];
         }
     }
 
@@ -834,16 +826,19 @@ int main(int argc, char *argv[])
                          << ", Lambda = " << eigenvalues[i] << endl;
                 }
 
-                Vector ev;
-                if (fom || offline) {
-                    ev = lobpcg->GetEigenvector(i);
+                Vector eigenvector_i;
+                if (fom || offline)
+                {
+                    eigenvector_i = lobpcg->GetEigenvector(i);
                 }
-                else {
-                    evect.GetRow(i, ev);
+                else
+                {
+                    eigenvectors.GetRow(i, eigenvector_i);
                 }
                 // convert eigenvector from HypreParVector to ParGridFunction
-                ev *= sign_ev[i] / sqrt(InnerProduct(ev, ev));
-                x = ev;
+                eigenvector_i *= sign_eigenvectors[i] / sqrt(InnerProduct(eigenvector_i,
+                                 eigenvector_i));
+                x = eigenvector_i;
 
                 sout << "parallel " << num_procs << " " << myid << "\n"
                      << "solution\n" << *pmesh << x << flush
@@ -910,14 +905,18 @@ int main(int argc, char *argv[])
 
 double Conductivity(const Vector &x)
 {
-    Vector center(x.Size());
     switch (problem)
     {
     case 1:
         return 1.0;
     case 2:
-        return 1.0 + gaussian_depth * std::exp(-x.DistanceSquaredTo(center) / pow(
-                gaussian_width, 2.0));
+        double cx = 1.0 + amplitude;
+        for (int i = 0; i < x.Size(); ++i)
+        {
+            if (8 * abs(x(i) - center(i)) > (bb_max[i] - bb_min[i]))
+                cx = 1.0;
+        }
+        return cx;
     case 3:
     case 4:
         return 1.0;
@@ -927,16 +926,17 @@ double Conductivity(const Vector &x)
 
 double Potential(const Vector &x)
 {
+    double d_sq = x.DistanceSquaredTo(center);
     switch (problem)
     {
     case 1:
     case 2:
         return 0.0;
     case 3:
-        return gaussian_depth * std::exp(-x.DistanceSquaredTo(center) / pow(
-                                             gaussian_width, 2.0));
+        double radius = 5.0 * h_max;
+        return amplitude * std::exp(-d_sq / pow(radius, 2.0));
     case 4:
-        return gaussian_depth / x.DistanceSquaredTo(center);
+        return amplitude / d_sq;
     }
     return 0.0;
 }

From 88a9368d51ba25ab5591864fa5af1cada4715483 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 20 Jun 2024 07:50:10 -0700
Subject: [PATCH 26/50] Remove verbose

---
 .../prom/elliptic_eigenproblem_global_rom.cpp    | 16 ----------------
 1 file changed, 16 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index ab32723e6..086b88b04 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -729,25 +729,9 @@ int main(int argc, char *argv[])
                     mode_rom_ifs.precision(16);
                     mode_rom.Load(mode_rom_ifs, eigenvectors.NumCols());
                     mode_rom_ifs.close();
-                    std::cout << "Projecting FOM mode " << i
-                              << " onto ROM mode " << j << std::endl;
-                    std::cout << "ip = " <<  InnerProduct(mode_fom, mode_rom) << std::endl;
                     mode_fom.Add(-InnerProduct(mode_fom, mode_rom), mode_rom);
                     mode_name.str("");
                 }
-                else
-                {
-                    mode_name << "mode_rom_" << setfill('0') << setw(2) << j << "."
-                              << setfill('0') << setw(6) << myid;
-                    ifstream mode_rom_ifs(mode_name.str().c_str());
-                    mode_rom_ifs.precision(16);
-                    mode_rom.Load(mode_rom_ifs, eigenvectors.NumCols());
-                    mode_rom_ifs.close();
-                    std::cout << "Not projecting FOM mode " << i
-                              << " onto ROM mode " << j << std::endl;
-                    std::cout << "ip = " <<  InnerProduct(mode_fom, mode_rom) << std::endl;
-                    mode_name.str("");
-                }
             }
 
             const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom));

From 60f0013c21689f673dbf46a8e485a2c61674604f Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 20 Jun 2024 08:53:05 -0700
Subject: [PATCH 27/50] Update example in header

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 54 +++++++++----------
 1 file changed, 26 insertions(+), 28 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 086b88b04..2c7810f07 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -23,39 +23,37 @@
 //               order < 1 (quadratic for quadratic curvilinear mesh, NURBS for
 //               NURBS mesh, etc.)
 //
-// Offline phase: elliptic_eigenproblem_global_rom -offline -p 2 -id 0 -a 0.40
-//                elliptic_eigenproblem_global_rom -offline -p 2 -id 1 -a 0.45
-//                elliptic_eigenproblem_global_rom -offline -p 2 -id 2 -a 0.55
-//                elliptic_eigenproblem_global_rom -offline -p 2 -id 3 -a 0.60
+// Offline phase: elliptic_eigenproblem_global_rom -offline -p 2 -rs 2 -n 4 -id 0 -a 0
+//                elliptic_eigenproblem_global_rom -offline -p 2 -rs 2 -n 4 -id 1 -a 1
 //
-// Merge phase:   elliptic_eigenproblem_global_rom -merge -p 2 -ns 4
+// Merge phase:   elliptic_eigenproblem_global_rom -merge -p 2 -rs 2 -n 4 -ns 2
 //
 // FOM run (for error calculation):
-//                elliptic_eigenproblem_global_rom -fom -p 2 -a 0.50
-// Example output:
-//   Number of unknowns: 289
-//    Eigenvalue 0: 0.04533314
-//    Eigenvalue 1: 0.11332411
-//    Eigenvalue 2: 0.11486387
-//    Eigenvalue 3: 0.18192
-//    Eigenvalue 4: 0.22964377
-//   Elapsed time for assembling FOM: 1.471708e-03 second
-//   Elapsed time for solving FOM: 3.416310e-01 second
+//                elliptic_eigenproblem_global_rom -fom -p 2 -rs 2 -n 4 -a 0.5
+//
+// Online phase:  elliptic_eigenproblem_global_rom -online -p 2 -rs 2 -n 4 -a 0.5 -ef 1.0
 //
-// Online phase:  elliptic_eigenproblem_global_rom -online -p 2 -a 0.50
 // Example output:
-//   Eigenvalue 0: = 0.048430949
-//   Eigenvalue 1: = 0.12021157
-//   Eigenvalue 2: = 0.12147847
-//   Eigenvalue 3: = 0.19456504
-//   Eigenvalue 4: = 0.24285855
-//   Relative error of ROM solution for eigenvalue 0 = 0.068334316
-//   Relative error of ROM solution for eigenvalue 1 = 0.060776586
-//   Relative error of ROM solution for eigenvalue 2 = 0.057586388
-//   Relative error of ROM solution for eigenvalue 3 = 0.069508795
-//   Relative error of ROM solution for eigenvalue 4 = 0.057544708
-//   Elapsed time for assembling ROM: 4.289041e-03 second
-//   Elapsed time for solving ROM: 6.225410e-04 second
+//   FOM solution for eigenvalue 0 = 19.878564
+//   ROM solution for eigenvalue 0 = 19.878631
+//   Absolute error of ROM solution for eigenvalue 0 = 6.6618756e-05
+//   Relative error of ROM solution for eigenvalue 0 = 3.3512861e-06
+//   FOM solution for eigenvalue 1 = 53.177886
+//   ROM solution for eigenvalue 1 = 53.179662
+//   Absolute error of ROM solution for eigenvalue 1 = 0.0017768914
+//   Relative error of ROM solution for eigenvalue 1 = 3.3414104e-05
+//   FOM solution for eigenvalue 2 = 53.177886
+//   ROM solution for eigenvalue 2 = 53.179662
+//   Absolute error of ROM solution for eigenvalue 2 = 0.0017769039
+//   Relative error of ROM solution for eigenvalue 2 = 3.3414339e-05
+//   FOM solution for eigenvalue 3 = 81.245811
+//   ROM solution for eigenvalue 3 = 81.246105
+//   Absolute error of ROM solution for eigenvalue 3 = 0.00029387266
+//   Relative error of ROM solution for eigenvalue 3 = 3.6170808e-06
+//   Relative l2 error of ROM eigenvector 0 = 0.0023815683
+//   Relative l2 error of ROM eigenvector 1 = 0.0097810599
+//   Relative l2 error of ROM eigenvector 2 = 0.0097816076
+//   Relative l2 error of ROM eigenvector 3 = 0.0043291658
 
 #include "mfem.hpp"
 #include "linalg/BasisGenerator.h"

From 4d0abab406a5ebc685c6b6c6bf774ca4df3c4140 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 20 Jun 2024 09:02:38 -0700
Subject: [PATCH 28/50] Fix failing CI

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 2c7810f07..73b6388b0 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -887,12 +887,13 @@ int main(int argc, char *argv[])
 
 double Conductivity(const Vector &x)
 {
+    double cx;
     switch (problem)
     {
     case 1:
         return 1.0;
     case 2:
-        double cx = 1.0 + amplitude;
+        cx = 1.0 + amplitude;
         for (int i = 0; i < x.Size(); ++i)
         {
             if (8 * abs(x(i) - center(i)) > (bb_max[i] - bb_min[i]))
@@ -908,6 +909,7 @@ double Conductivity(const Vector &x)
 
 double Potential(const Vector &x)
 {
+    double radius = 5.0 * h_max;
     double d_sq = x.DistanceSquaredTo(center);
     switch (problem)
     {
@@ -915,10 +917,9 @@ double Potential(const Vector &x)
     case 2:
         return 0.0;
     case 3:
-        double radius = 5.0 * h_max;
         return amplitude * std::exp(-d_sq / pow(radius, 2.0));
     case 4:
-        return amplitude / d_sq;
+        return amplitude * d_sq;
     }
     return 0.0;
 }

From 16fc913d994cadb5c25220ed0c74408e032a7d2c Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Tue, 9 Jul 2024 09:06:12 -0700
Subject: [PATCH 29/50] Updates

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 378 ++++++++----------
 1 file changed, 171 insertions(+), 207 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 73b6388b0..3d55f5ee9 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -71,8 +71,8 @@ double Conductivity(const Vector &x);
 double Potential(const Vector &x);
 int problem = 1;
 double amplitude = -800.0;
-double relative_center = 0.0;
-Vector center;
+double pseudo_time = 0.0;
+int potential_well_switch = 0;
 Vector bb_min, bb_max;
 double h_min, h_max, k_min, k_max;
 
@@ -105,7 +105,7 @@ int main(int argc, char *argv[])
     bool online = false;
     int id = 0;
     int nsets = 0;
-    double ef = 0.9999;
+    double ef = 1.0;
     int rdim = -1;
     int verbose_level = 0;
 
@@ -131,8 +131,11 @@ int main(int argc, char *argv[])
                    "Random seed used to initialize LOBPCG.");
     args.AddOption(&amplitude, "-a", "--amplitude",
                    "Amplitude of coefficient fields.");
-    args.AddOption(&relative_center, "-c", "--center",
-                   "Number of grid elements to which center is shifted.");
+    args.AddOption(&pseudo_time, "-t", "--pseudo-time",
+                   "Pseudo-time of the motion of the center.");
+    args.AddOption(&potential_well_switch, "-pw", "--potential-well",
+                   "Customized which potential well is turned off in problem 4.\n\t"
+                   "Available options: 0 (default, no wells off), 1, 2.");
     args.AddOption(&id, "-id", "--id", "Parametric id");
     args.AddOption(&nsets, "-ns", "--nset", "Number of parametric snapshot sets");
     args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -227,12 +230,6 @@ int main(int argc, char *argv[])
     mesh->GetCharacteristics(h_min, h_max, k_min, k_max);
     h_max *= pow(0.5, ser_ref_levels + par_ref_levels);
 
-    center.SetSize(dim);
-    for (int i = 0; i < dim; i++)
-    {
-        center(i) = h_max * relative_center + 0.5 * (bb_min[i] + bb_max[i]);
-    }
-
     // 4. Refine the mesh in serial to increase the resolution. In this example
     //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
     //    a command-line parameter.
@@ -334,7 +331,7 @@ int main(int argc, char *argv[])
     FunctionCoefficient kappa_0(Conductivity);
     FunctionCoefficient v_0(Potential);
 
-    // Project conductivity and potential to visualize initial condition
+    // Project conductivity and potential for visualization
     ParGridFunction c_gf(fespace);
     c_gf.ProjectCoefficient(kappa_0);
 
@@ -386,13 +383,14 @@ int main(int argc, char *argv[])
     delete a;
     delete m;
 
-    ParGridFunction x(fespace);
+    ParGridFunction eigenfunction_i(fespace);
+    Vector eigenvector_i(size);
     HypreLOBPCG *lobpcg;
     Array<double> eigenvalues;
-    DenseMatrix eigenvectors;
+    std::vector<Vector> eigenvectors;
 
     // 11. The offline phase
-    if(fom || offline)
+    if (fom || offline)
     {
         // 12. Define and configure the LOBPCG eigensolver and the BoomerAMG
         //     preconditioner for A to be used within the solver. Set the matrices
@@ -472,8 +470,9 @@ int main(int argc, char *argv[])
             }
             if (offline)
             {
-                x = lobpcg->GetEigenvector(i);
-                generator->takeSample(x.GetData());
+                eigenfunction_i = lobpcg->GetEigenvector(i);
+                eigenfunction_i /= sqrt(InnerProduct(eigenfunction_i, eigenfunction_i));
+                generator->takeSample(eigenfunction_i.GetData());
             }
         }
 
@@ -494,7 +493,6 @@ int main(int argc, char *argv[])
         assembleTimer.Start();
         CAROM::BasisReader reader(basis_filename);
 
-        Vector eigenvalues_rom;
         const CAROM::Matrix *spatialbasis;
         if (rdim != -1)
         {
@@ -509,6 +507,7 @@ int main(int argc, char *argv[])
         const int numColumnRB = spatialbasis->numColumns();
         if (myid == 0) printf("spatial basis dimension is %d x %d\n", numRowRB,
                                   numColumnRB);
+        MFEM_VERIFY(numColumnRB >= nev, "basis dimension >= number of eigenvectors");
 
         // 17. form ROM operator
         CAROM::Matrix ReducedA;
@@ -525,219 +524,166 @@ int main(int argc, char *argv[])
 
         assembleTimer.Stop();
 
+        Vector eigenvalues_rom;
+        DenseMatrix reduced_eigenvectors;
         // 18. solve ROM
         solveTimer.Start();
         // (Q^T A Q) c = \lambda (Q^T M Q) c
-        A_mat->Eigenvalues(*M_mat, eigenvalues_rom, eigenvectors);
+        A_mat->Eigenvalues(*M_mat, eigenvalues_rom, reduced_eigenvectors);
         solveTimer.Stop();
 
-        if (myid == 0)
+        // Postprocessing
+        eigenvalues = Array<double>(eigenvalues_rom.GetData(), eigenvalues_rom.Size());
+        Vector mode_rom(numRowRB);
+        std::vector<Vector> modes_rom;
+        for (int j = 0; j < nev; j++)
+        {
+            Vector reduced_eigenvector_j;
+            reduced_eigenvectors.GetColumn(j, reduced_eigenvector_j);
+            CAROM::Vector reduced_eigenvector_j_carom(reduced_eigenvector_j.GetData(),
+                    reduced_eigenvector_j.Size(), false, false);
+            CAROM::Vector *eigenvector_j_carom = spatialbasis->mult(
+                    reduced_eigenvector_j_carom);
+            mode_rom.SetData(eigenvector_j_carom->getData());
+            mode_rom /= sqrt(InnerProduct(mode_rom, mode_rom));
+            modes_rom.push_back(mode_rom);
+            delete eigenvector_j_carom;
+        }
+
+        // Calculate errors of eigenvalues
+        ostringstream sol_eigenvalue_name_fom;
+        sol_eigenvalue_name_fom << "sol_eigenvalues_fom." << setfill('0') << setw(
+                                    6) << myid;
+        Vector eigenvalues_fom(nev);
+        ifstream fom_file;
+        fom_file.open(sol_eigenvalue_name_fom.str().c_str());
+        eigenvalues_fom.Load(fom_file, eigenvalues_fom.Size());
+        fom_file.close();
+
+        Vector diff_eigenvalues(nev);
+        for (int i = 0; i < nev; i++)
         {
-            eigenvalues = Array<double>(eigenvalues_rom.GetData(), eigenvalues_rom.Size());
-            for (int i = 0; i < eigenvalues_rom.Size() && i < nev; i++)
+            diff_eigenvalues[i] = eigenvalues_fom[i] - eigenvalues_rom[i];
+            if (myid == 0)
             {
-                std::cout << "Eigenvalue " << i << ": = " << eigenvalues[i] << "\n";
+                std::cout << "FOM solution for eigenvalue " << i << " = " <<
+                          eigenvalues_fom[i] << std::endl;
+                std::cout << "ROM solution for eigenvalue " << i << " = " <<
+                          eigenvalues_rom[i] << std::endl;
+                std::cout << "Absolute error of ROM solution for eigenvalue " << i << " = " <<
+                          abs(diff_eigenvalues[i]) << std::endl;
+                std::cout << "Relative error of ROM solution for eigenvalue " << i << " = " <<
+                          abs(diff_eigenvalues[i]) / abs(eigenvalues_fom[i]) << std::endl;
             }
         }
 
-        DenseMatrix tmp(eigenvectors);
-        eigenvectors = DenseMatrix(nev, numRowRB);
-        for (int i = 0; i < eigenvalues_rom.Size() && i < nev; i++)
+        // Calculate errors of eigenvectors
+        ostringstream mode_name_fom;
+        Vector mode_fom(numRowRB);
+        for (int i = 0; i < nev; i++)
         {
-            Vector reduced_eigenvector;
-            tmp.GetRow(i, reduced_eigenvector);
-            CAROM::Vector reduced_eigenvector_carom(reduced_eigenvector.GetData(),
-                                                    reduced_eigenvector.Size(), false, false);
-            CAROM::Vector *eigenvector_carom = spatialbasis->mult(
-                                                   reduced_eigenvector_carom);
-            eigenvectors.SetRow(i, eigenvector_carom->getData());
-            delete eigenvector_carom;
+            mode_name_fom << "eigenfunction_fom_" << setfill('0') << setw(2) << i << "."
+                          << setfill('0') << setw(6) << myid;
+
+            ifstream mode_fom_ifs(mode_name_fom.str().c_str());
+            mode_fom_ifs.precision(16);
+            mode_fom.Load(mode_fom_ifs, numRowRB);
+            mode_fom_ifs.close();
+
+            eigenvector_i = 0.0;
+            for (int j = 0; j < nev; j++)
+            {
+                if (myid == 0)
+                {
+                    std::cout << "correlation_matrix(" << i+1 << "," << j+1 << ") = " <<
+                              InnerProduct(mode_fom, modes_rom[j]) << std::endl;
+                }
+                if (abs(eigenvalues_fom[j] - eigenvalues_fom[i]) < 1e-6)
+                {
+                    eigenvector_i.Add(InnerProduct(mode_fom, modes_rom[j]), modes_rom[j]);
+                }
+            }
+            //eigenvector_i /= sqrt(InnerProduct(eigenvector_i, eigenvector_i));
+            eigenvectors.push_back(eigenvector_i);
+            mode_fom.Add(-1.0, eigenvector_i);
+            const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom));
+            if (myid == 0) std::cout << "Relative l2 error of ROM eigenvector " << i <<
+                                         " = " << diffNorm << std::endl;
+
+            mode_name_fom.str("");
         }
 
         delete spatialbasis;
-
         delete A_mat;
         delete M_mat;
     }
 
-    ostringstream sol_eigenvalue_name, sol_eigenvalue_name_fom;
-    if (fom || offline)
-    {
-        sol_eigenvalue_name << "sol_eigenvalues_fom." << setfill('0') << setw(
-                                6) << myid;
-    }
-    if (online)
-    {
-        sol_eigenvalue_name << "sol_eigenvalues." << setfill('0') << setw(6) << myid;
-        sol_eigenvalue_name_fom << "sol_eigenvalues_fom." << setfill('0') << setw(
-                                    6) << myid;
-    }
-
     // 19. Save the refined mesh and the modes in parallel. This output can be
     //     viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
-    Vector sign_eigenvectors(nev);
     {
-        ostringstream mesh_name, mode_name, mode_ref_name;
+        ostringstream sol_eigenvalue_name;
+        if (fom || offline)
+        {
+            sol_eigenvalue_name << "sol_eigenvalues_fom." << setfill('0') << setw(
+                                    6) << myid;
+        }
+        if (online)
+        {
+            sol_eigenvalue_name << "sol_eigenvalues." << setfill('0') << setw(6) << myid;
+        }
+
+        ofstream sol_eigenvalue_ofs(sol_eigenvalue_name.str().c_str());
+        sol_eigenvalue_ofs.precision(16);
+        for (int i = 0; i < nev; ++i)
+        {
+            sol_eigenvalue_ofs << eigenvalues[i] << std::endl;
+        }
+
+        ostringstream mesh_name, mode_name;
         mesh_name << "elliptic_eigenproblem-mesh." << setfill('0') << setw(6) << myid;
 
         ofstream mesh_ofs(mesh_name.str().c_str());
         mesh_ofs.precision(8);
         pmesh->Print(mesh_ofs);
 
-        std::string mode_prefix = "mode_";
-        std::string mode_ref_prefix = "mode_";
+        std::string mode_prefix = "eigenfunction_";
         if (fom || offline)
         {
             mode_prefix += "fom_";
-            mode_ref_prefix += "ref_";
         }
         else if (online)
         {
             mode_prefix += "rom_";
-            mode_ref_prefix += "fom_";
         }
 
-        for (int i=0; i < nev && i < eigenvalues.Size(); i++)
+        for (int i = 0; i < nev; i++)
         {
-            Vector eigenvector_i;
             if (fom || offline)
             {
-                eigenvector_i = lobpcg->GetEigenvector(i);
+                eigenfunction_i = lobpcg->GetEigenvector(i);
+                eigenfunction_i /= sqrt(InnerProduct(eigenfunction_i, eigenfunction_i));
             }
             else
             {
                 // for online, eigenvectors are stored in eigenvectors matrix
-                eigenvectors.GetRow(i, eigenvector_i);
-            }
-
-            Vector mode_ref(eigenvector_i.Size());
-            mode_ref_name << mode_ref_prefix << setfill('0') << setw(2) << i << "."
-                          << setfill('0') << setw(6) << myid;
-            if (offline && id == 0)
-            {
-                mode_ref = eigenvector_i;
-                ofstream mode_ref_ofs(mode_ref_name.str().c_str());
-                mode_ref_ofs.precision(16);
-
-                // TODO: issue using .Load() function if file written with .Save()?
-                //x.Save(mode_ref_ofs);
-                for (int j = 0; j < x.Size(); j++)
-                {
-                    mode_ref_ofs << mode_ref[j] << "\n";
-                }
-                mode_ref_ofs.close();
-            }
-            else
-            {
-                ifstream mode_ref_ifs(mode_ref_name.str().c_str());
-                mode_ref_ifs.precision(16);
-                mode_ref.Load(mode_ref_ifs, eigenvector_i.Size());
-                mode_ref_ifs.close();
-            }
-            mode_ref_name.str("");
-            sign_eigenvectors[i] = (InnerProduct(mode_ref, eigenvector_i) >= 0) ? 1 : -1;
-            eigenvector_i *= sign_eigenvectors[i] / sqrt(InnerProduct(eigenvector_i,
-                             eigenvector_i));
-
-            if (InnerProduct(mode_ref, eigenvector_i) < 0.9)
-            {
-                std::cout << "Warning: eigenvector " << i <<
-                          " in FOM and ROM are not directly comparable."
-                          << std::endl;
-                std::cout << "Inner product = " << InnerProduct(mode_ref,
-                          eigenvector_i) << std::endl;
-                std::cout << "TODO: Visualization of projected eigenvector." << std::endl;
+                // convert eigenvector from HypreParVector to ParGridFunction
+                eigenfunction_i = eigenvectors[i];
             }
 
-            // convert eigenvector from HypreParVector to ParGridFunction
-            x = eigenvector_i;
             mode_name << mode_prefix << setfill('0') << setw(2) << i << "."
                       << setfill('0') << setw(6) << myid;
 
             ofstream mode_ofs(mode_name.str().c_str());
             mode_ofs.precision(16);
             // TODO: issue using .Load() function if file written with .Save()?
-            //x.Save(mode_ofs);
-            for (int j = 0; j < x.Size(); j++)
+            //eigenfunction_i.Save(mode_ofs);
+            for (int j = 0; j < eigenfunction_i.Size(); j++)
             {
-                mode_ofs << x[j] << "\n";
+                mode_ofs << eigenfunction_i[j] << "\n";
             }
             mode_ofs.close();
             mode_name.str("");
         }
-
-        ofstream sol_eigenvalue_ofs(sol_eigenvalue_name.str().c_str());
-        sol_eigenvalue_ofs.precision(16);
-        for (int i = 0; i < nev && i < eigenvalues.Size(); ++i)
-        {
-            sol_eigenvalue_ofs << eigenvalues[i] << std::endl;
-        }
-    }
-
-    if (online)
-    {
-        // Initialize FOM solution
-        Vector eigenvalues_fom(nev);
-
-        ifstream fom_file;
-        fom_file.open(sol_eigenvalue_name_fom.str().c_str());
-        eigenvalues_fom.Load(fom_file, eigenvalues_fom.Size());
-        fom_file.close();
-
-        Vector diff_eigenvalues(nev);
-        for (int i = 0; i < eigenvalues.Size() && i < nev; i++)
-        {
-            diff_eigenvalues[i] = eigenvalues_fom[i] - eigenvalues[i];
-            if (myid == 0)
-            {
-                std::cout << "FOM solution for eigenvalue " << i << " = " <<
-                          eigenvalues_fom[i] << std::endl;
-                std::cout << "ROM solution for eigenvalue " << i << " = " <<
-                          eigenvalues[i] << std::endl;
-                std::cout << "Absolute error of ROM solution for eigenvalue " << i << " = " <<
-                          abs(diff_eigenvalues[i]) << std::endl;
-                std::cout << "Relative error of ROM solution for eigenvalue " << i << " = " <<
-                          abs(diff_eigenvalues[i]) / abs(eigenvalues_fom[i]) << std::endl;
-            }
-        }
-
-        // Calculate errors of eigenvectors
-        ostringstream mode_name, mode_name_fom;
-        Vector mode_rom(eigenvectors.NumCols());
-        Vector mode_fom(eigenvectors.NumCols());
-        for (int i = 0; i < eigenvalues.Size() && i < nev; i++)
-        {
-            mode_name_fom << "mode_fom_" << setfill('0') << setw(2) << i << "."
-                          << setfill('0') << setw(6) << myid;
-
-            ifstream mode_fom_ifs(mode_name_fom.str().c_str());
-            mode_fom_ifs.precision(16);
-            mode_fom.Load(mode_fom_ifs, eigenvectors.NumCols());
-            mode_fom_ifs.close();
-
-            const double fomNorm = sqrt(InnerProduct(mode_fom, mode_fom));
-
-            for (int j = 0; j < eigenvalues.Size() && j < nev; j++)
-            {
-                if (abs(eigenvalues_fom[j] - eigenvalues_fom[i]) < 1e-6)
-                {
-                    mode_name << "mode_rom_" << setfill('0') << setw(2) << j << "."
-                              << setfill('0') << setw(6) << myid;
-                    ifstream mode_rom_ifs(mode_name.str().c_str());
-                    mode_rom_ifs.precision(16);
-                    mode_rom.Load(mode_rom_ifs, eigenvectors.NumCols());
-                    mode_rom_ifs.close();
-                    mode_fom.Add(-InnerProduct(mode_fom, mode_rom), mode_rom);
-                    mode_name.str("");
-                }
-            }
-
-            const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom));
-            if (myid == 0) std::cout << "Relative l2 error of ROM eigenvector " << i <<
-                                         " = " << diffNorm / fomNorm << std::endl;
-
-            mode_name_fom.str("");
-        }
     }
 
     VisItDataCollection *visit_dc = NULL;
@@ -750,22 +696,21 @@ int main(int argc, char *argv[])
         visit_dc->RegisterField("Conductivity", &c_gf);
         visit_dc->RegisterField("Potential", &p_gf);
         std::vector<ParGridFunction*> visit_eigenvectors;
-        for (int i = 0; i < nev && i < eigenvalues.Size(); i++)
+        for (int i = 0; i < nev; i++)
         {
-            Vector eigenvector_i;
             if (fom || offline)
             {
-                eigenvector_i = lobpcg->GetEigenvector(i);
+                eigenfunction_i = lobpcg->GetEigenvector(i);
+                eigenfunction_i /= sqrt(InnerProduct(eigenfunction_i, eigenfunction_i));
             }
-            else {
+            else
+            {
                 // for online, eigenvectors are stored in eigvenvectors matrix
-                eigenvectors.GetRow(i, eigenvector_i);
+                // convert eigenvector from HypreParVector to ParGridFunction
+                eigenfunction_i = eigenvectors[i];
             }
-            eigenvector_i *= sign_eigenvectors[i] / sqrt(InnerProduct(eigenvector_i,
-                             eigenvector_i));
-            // convert eigenvector from HypreParVector to ParGridFunction
-            x = eigenvector_i;
-            visit_eigenvectors.push_back(new ParGridFunction(x));
+
+            visit_eigenvectors.push_back(new ParGridFunction(eigenfunction_i));
             visit_dc->RegisterField("Eigenmode_" + std::to_string(i),
                                     visit_eigenvectors.back());
         }
@@ -800,7 +745,7 @@ int main(int argc, char *argv[])
         }
         else
         {
-            for (int i = 0; i < nev && i < eigenvalues.Size(); i++)
+            for (int i = 0; i < nev; i++)
             {
                 if ( myid == 0 )
                 {
@@ -808,22 +753,20 @@ int main(int argc, char *argv[])
                          << ", Lambda = " << eigenvalues[i] << endl;
                 }
 
-                Vector eigenvector_i;
                 if (fom || offline)
                 {
-                    eigenvector_i = lobpcg->GetEigenvector(i);
+                    eigenfunction_i = lobpcg->GetEigenvector(i);
+                    eigenfunction_i /= sqrt(InnerProduct(eigenfunction_i, eigenfunction_i));
                 }
                 else
                 {
-                    eigenvectors.GetRow(i, eigenvector_i);
+                    // for online, eigenvectors are stored in eigvenvectors matrix
+                    // convert eigenvector from HypreParVector to ParGridFunction
+                    eigenfunction_i = eigenvectors[i];
                 }
-                // convert eigenvector from HypreParVector to ParGridFunction
-                eigenvector_i *= sign_eigenvectors[i] / sqrt(InnerProduct(eigenvector_i,
-                                 eigenvector_i));
-                x = eigenvector_i;
 
                 sout << "parallel " << num_procs << " " << myid << "\n"
-                     << "solution\n" << *pmesh << x << flush
+                     << "solution\n" << *pmesh << eigenfunction_i << flush
                      << "window_title 'Eigenmode " << i+1 << '/' << nev
                      << ", Lambda = " << eigenvalues[i] << "'" << endl;
 
@@ -847,13 +790,13 @@ int main(int argc, char *argv[])
     // 20. print timing info
     if (myid == 0)
     {
-        if(fom || offline)
+        if (fom || offline)
         {
             printf("Elapsed time for assembling FOM: %e second\n",
                    assembleTimer.RealTime());
             printf("Elapsed time for solving FOM: %e second\n", solveTimer.RealTime());
         }
-        if(online)
+        if (online)
         {
             printf("Elapsed time for assembling ROM: %e second\n",
                    assembleTimer.RealTime());
@@ -887,19 +830,15 @@ int main(int argc, char *argv[])
 
 double Conductivity(const Vector &x)
 {
-    double cx;
     switch (problem)
     {
     case 1:
         return 1.0;
     case 2:
-        cx = 1.0 + amplitude;
         for (int i = 0; i < x.Size(); ++i)
-        {
-            if (8 * abs(x(i) - center(i)) > (bb_max[i] - bb_min[i]))
-                cx = 1.0;
-        }
-        return cx;
+            if (8 * abs(x(i) - 0.5 * (bb_max[i] + bb_min[i])) > bb_max[i] - bb_min[i])
+                return 1.0;
+        return 1.0 + amplitude;
     case 3:
     case 4:
         return 1.0;
@@ -909,17 +848,42 @@ double Conductivity(const Vector &x)
 
 double Potential(const Vector &x)
 {
-    double radius = 5.0 * h_max;
-    double d_sq = x.DistanceSquaredTo(center);
+    Vector center(x.Size());
+    double rho = 0.0;
     switch (problem)
     {
     case 1:
     case 2:
         return 0.0;
     case 3:
-        return amplitude * std::exp(-d_sq / pow(radius, 2.0));
+        // center = (0.5 + t*h, 0.5)
+        center(0) = 0.5 * (bb_min[0] + bb_max[0]) + pseudo_time * h_max;
+        for (int i = 1; i < x.Size(); i++)
+            center(i) = 0.5 * (bb_min[i] + bb_max[i]);
+        return amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,
+                                    2.0));
     case 4:
-        return amplitude * d_sq;
+        if (potential_well_switch != 1)
+        {
+            // add well with center = (1/4 + cos(2*pi*t) / 64, 1/4 + sin(2*pi*t) / 64, 1/4)
+            center(0) = 0.25 * (bb_min[0] + bb_max[0]) + cos(2 * M_PI * pseudo_time) / 64;
+            center(1) = 0.25 * (bb_min[1] + bb_max[1]) + sin(2 * M_PI * pseudo_time) / 64;
+            for (int i = 2; i < x.Size(); i++)
+                center(i) = 0.25 * (bb_min[i] + bb_max[i]);
+            rho += amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,
+                                        2.0));
+        }
+        if (potential_well_switch != 2)
+        {
+            // add well with center = (3/4 + sin(2*pi*t) / 64, 3/4 + cos(2*pi*t) / 64, 3/4)
+            center(0) = 0.75 * (bb_min[0] + bb_max[0]) - cos(2 * M_PI * pseudo_time) / 64;
+            center(1) = 0.75 * (bb_min[1] + bb_max[1]) - sin(2 * M_PI * pseudo_time) / 64;
+            for (int i = 2; i < x.Size(); i++)
+                center(i) = 0.75 * (bb_min[i] + bb_max[i]);
+            rho += amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,
+                                        2.0));
+        }
+        return rho;
     }
     return 0.0;
 }

From a1999253565e950af1f2747e0e32d3e9eb42e15d Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Tue, 9 Jul 2024 14:55:17 -0700
Subject: [PATCH 30/50] Fix verbose

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 3d55f5ee9..955580ff9 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -593,10 +593,10 @@ int main(int argc, char *argv[])
             eigenvector_i = 0.0;
             for (int j = 0; j < nev; j++)
             {
-                if (myid == 0)
+                if (myid == 0 && verbose_level > 0)
                 {
                     std::cout << "correlation_matrix(" << i+1 << "," << j+1 << ") = " <<
-                              InnerProduct(mode_fom, modes_rom[j]) << std::endl;
+                              InnerProduct(mode_fom, modes_rom[j]) << ";" << std::endl;
                 }
                 if (abs(eigenvalues_fom[j] - eigenvalues_fom[i]) < 1e-6)
                 {
@@ -866,8 +866,8 @@ double Potential(const Vector &x)
         if (potential_well_switch != 1)
         {
             // add well with center = (1/4 + cos(2*pi*t) / 64, 1/4 + sin(2*pi*t) / 64, 1/4)
-            center(0) = 0.25 * (bb_min[0] + bb_max[0]) + cos(2 * M_PI * pseudo_time) / 64;
-            center(1) = 0.25 * (bb_min[1] + bb_max[1]) + sin(2 * M_PI * pseudo_time) / 64;
+            center(0) = (0.25 + cos(2 * M_PI * pseudo_time) / 64) * (bb_min[0] + bb_max[0]);
+            center(1) = (0.25 + sin(2 * M_PI * pseudo_time) / 64) * (bb_min[1] + bb_max[1]);
             for (int i = 2; i < x.Size(); i++)
                 center(i) = 0.25 * (bb_min[i] + bb_max[i]);
             rho += amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,
@@ -876,8 +876,8 @@ double Potential(const Vector &x)
         if (potential_well_switch != 2)
         {
             // add well with center = (3/4 + sin(2*pi*t) / 64, 3/4 + cos(2*pi*t) / 64, 3/4)
-            center(0) = 0.75 * (bb_min[0] + bb_max[0]) - cos(2 * M_PI * pseudo_time) / 64;
-            center(1) = 0.75 * (bb_min[1] + bb_max[1]) - sin(2 * M_PI * pseudo_time) / 64;
+            center(0) = (0.75 - cos(2 * M_PI * pseudo_time) / 64) * (bb_min[0] + bb_max[0]);
+            center(1) = (0.75 - sin(2 * M_PI * pseudo_time) / 64) * (bb_min[1] + bb_max[1]);
             for (int i = 2; i < x.Size(); i++)
                 center(i) = 0.75 * (bb_min[i] + bb_max[i]);
             rho += amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,

From 2505aaa14a1c142be25c52fd8f9b9264f34b43a1 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Tue, 9 Jul 2024 19:33:23 -0700
Subject: [PATCH 31/50] Minor changes to verbose

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 19 +++++++++----------
 1 file changed, 9 insertions(+), 10 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 955580ff9..7f29f4064 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -595,7 +595,7 @@ int main(int argc, char *argv[])
             {
                 if (myid == 0 && verbose_level > 0)
                 {
-                    std::cout << "correlation_matrix(" << i+1 << "," << j+1 << ") = " <<
+                    std::cout << "correlation_matrix(" << j+1 << "," << i+1 << ") = " <<
                               InnerProduct(mode_fom, modes_rom[j]) << ";" << std::endl;
                 }
                 if (abs(eigenvalues_fom[j] - eigenvalues_fom[i]) < 1e-6)
@@ -603,7 +603,6 @@ int main(int argc, char *argv[])
                     eigenvector_i.Add(InnerProduct(mode_fom, modes_rom[j]), modes_rom[j]);
                 }
             }
-            //eigenvector_i /= sqrt(InnerProduct(eigenvector_i, eigenvector_i));
             eigenvectors.push_back(eigenvector_i);
             mode_fom.Add(-1.0, eigenvector_i);
             const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom));
@@ -863,21 +862,21 @@ double Potential(const Vector &x)
         return amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,
                                     2.0));
     case 4:
-        if (potential_well_switch != 1)
+        if (potential_well_switch == 0 || potential_well_switch == 1)
         {
-            // add well with center = (1/4 + cos(2*pi*t) / 64, 1/4 + sin(2*pi*t) / 64, 1/4)
-            center(0) = (0.25 + cos(2 * M_PI * pseudo_time) / 64) * (bb_min[0] + bb_max[0]);
-            center(1) = (0.25 + sin(2 * M_PI * pseudo_time) / 64) * (bb_min[1] + bb_max[1]);
+            // add well with first center
+            center(0) = (0.25 + 0.02 * cos(2 * M_PI * pseudo_time)) * (bb_min[0] + bb_max[0]);
+            center(1) = (0.25 + 0.02 * sin(2 * M_PI * pseudo_time)) * (bb_min[1] + bb_max[1]);
             for (int i = 2; i < x.Size(); i++)
                 center(i) = 0.25 * (bb_min[i] + bb_max[i]);
             rho += amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,
                                         2.0));
         }
-        if (potential_well_switch != 2)
+        if (potential_well_switch == 0 || potential_well_switch == 2)
         {
-            // add well with center = (3/4 + sin(2*pi*t) / 64, 3/4 + cos(2*pi*t) / 64, 3/4)
-            center(0) = (0.75 - cos(2 * M_PI * pseudo_time) / 64) * (bb_min[0] + bb_max[0]);
-            center(1) = (0.75 - sin(2 * M_PI * pseudo_time) / 64) * (bb_min[1] + bb_max[1]);
+            // add well with second center
+            center(0) = (0.75 - 0.02 * cos(2 * M_PI * pseudo_time)) * (bb_min[0] + bb_max[0]);
+            center(1) = (0.75 - 0.02 * sin(2 * M_PI * pseudo_time)) * (bb_min[1] + bb_max[1]);
             for (int i = 2; i < x.Size(); i++)
                 center(i) = 0.75 * (bb_min[i] + bb_max[i]);
             rho += amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,

From 825b330f0a746d02a8048dcd1c3d7a8f9dd8afde Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 11 Jul 2024 08:14:55 -0700
Subject: [PATCH 32/50] Astyle

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index e515649b0..e378345f2 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -896,8 +896,10 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
-            center(0) = (0.25 + 0.02 * cos(2 * M_PI * pseudo_time)) * (bb_min[0] + bb_max[0]);
-            center(1) = (0.25 + 0.02 * sin(2 * M_PI * pseudo_time)) * (bb_min[1] + bb_max[1]);
+            center(0) = (0.25 + 0.02 * cos(2 * M_PI * pseudo_time)) *
+                        (bb_min[0] + bb_max[0]);
+            center(1) = (0.25 + 0.02 * sin(2 * M_PI * pseudo_time)) *
+                        (bb_min[1] + bb_max[1]);
             for (int i = 2; i < x.Size(); i++)
                 center(i) = 0.25 * (bb_min[i] + bb_max[i]);
             rho += amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,
@@ -906,8 +908,10 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
-            center(0) = (0.75 - 0.02 * cos(2 * M_PI * pseudo_time)) * (bb_min[0] + bb_max[0]);
-            center(1) = (0.75 - 0.02 * sin(2 * M_PI * pseudo_time)) * (bb_min[1] + bb_max[1]);
+            center(0) = (0.75 - 0.02 * cos(2 * M_PI * pseudo_time)) *
+                        (bb_min[0] + bb_max[0]);
+            center(1) = (0.75 - 0.02 * sin(2 * M_PI * pseudo_time)) *
+                        (bb_min[1] + bb_max[1]);
             for (int i = 2; i < x.Size(); i++)
                 center(i) = 0.75 * (bb_min[i] + bb_max[i]);
             rho += amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,

From 42d09b7e2081bf1df419d2c5cafe97858d9741fc Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 11 Jul 2024 08:16:05 -0700
Subject: [PATCH 33/50] Clean up

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 31 -------------------
 1 file changed, 31 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index e378345f2..a36aaf768 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -619,7 +619,6 @@ int main(int argc, char *argv[])
 
     // 19. Save the refined mesh and the modes in parallel. This output can be
     //     viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
-    Vector sign_eigenvectors(nev);
     {
         ostringstream sol_eigenvalue_name;
         if (fom || offline)
@@ -650,12 +649,6 @@ int main(int argc, char *argv[])
         if (fom || offline)
         {
             mode_prefix += "fom_";
-            mode_ref_prefix += "ref_";
-        }
-        else if (online)
-        {
-            mode_prefix += "rom_";
-            mode_ref_prefix += "fom_";
         }
         else if (online)
         {
@@ -675,23 +668,7 @@ int main(int argc, char *argv[])
                 // convert eigenvector from HypreParVector to ParGridFunction
                 eigenfunction_i = eigenvectors[i];
             }
-            mode_ref_name.str("");
-            sign_eigenvectors[i] = (InnerProduct(mode_ref, eigenvector_i) >= 0) ? 1 : -1;
-            eigenvector_i *= sign_eigenvectors[i] / sqrt(InnerProduct(eigenvector_i,
-                             eigenvector_i));
 
-            if (InnerProduct(mode_ref, eigenvector_i) < 0.9)
-            {
-                std::cout << "Warning: eigenvector " << i <<
-                          " in FOM and ROM are not directly comparable."
-                          << std::endl;
-                std::cout << "Inner product = " << InnerProduct(mode_ref,
-                          eigenvector_i) << std::endl;
-                std::cout << "TODO: Visualization of projected eigenvector." << std::endl;
-            }
-
-            // convert eigenvector from HypreParVector to ParGridFunction
-            x = eigenvector_i;
             mode_name << mode_prefix << setfill('0') << setw(2) << i << "."
                       << setfill('0') << setw(6) << myid;
 
@@ -786,14 +763,6 @@ int main(int argc, char *argv[])
                     // convert eigenvector from HypreParVector to ParGridFunction
                     eigenfunction_i = eigenvectors[i];
                 }
-                else
-                {
-                    eigenvectors.GetRow(i, eigenvector_i);
-                }
-                // convert eigenvector from HypreParVector to ParGridFunction
-                eigenvector_i *= sign_eigenvectors[i] / sqrt(InnerProduct(eigenvector_i,
-                                 eigenvector_i));
-                x = eigenvector_i;
 
                 sout << "parallel " << num_procs << " " << myid << "\n"
                      << "solution\n" << *pmesh << eigenfunction_i << flush

From b42ab46bdd59ae7e81241bb633e8133a6a271ad9 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 18 Jul 2024 09:11:32 -0700
Subject: [PATCH 34/50] Make muber of LOBPCG iterations an option

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index a36aaf768..eaf902112 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -92,6 +92,7 @@ int main(int argc, char *argv[])
     int order = 1;
     int nev = 4;
     int seed = 75;
+    int lobpcg_niter = 200;
     bool slu_solver = false;
     bool sp_solver = false;
     bool cpardiso_solver = false;
@@ -160,6 +161,8 @@ int main(int argc, char *argv[])
                    "Reduced dimension.");
     args.AddOption(&verbose_level, "-v", "--verbose",
                    "Set the verbosity level of the LOBPCG solver and preconditioner. 0 is off.");
+    args.AddOption(&lobpcg_niter, "-fi", "--fom-iter",
+                   "Number of iterations for the LOBPCG solver.");
 #ifdef MFEM_USE_SUPERLU
     args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
                    "--no-superlu", "Use the SuperLU Solver.");

From 4742344fc405fcd7744f1ee4f33dbb3f0aa8f80c Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 8 Aug 2024 08:46:56 -0700
Subject: [PATCH 35/50] Set up the Gaussians for the Li2 example

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 27 ++++++++-----------
 1 file changed, 11 insertions(+), 16 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index eaf902112..95c0f5562 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -852,6 +852,7 @@ double Potential(const Vector &x)
 {
     Vector center(x.Size());
     double rho = 0.0;
+    double scaling_factor = pow((bb_max[0] - bb_min[0]) / 18.0, 2.0);
     switch (problem)
     {
     case 1:
@@ -868,26 +869,20 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
-            center(0) = (0.25 + 0.02 * cos(2 * M_PI * pseudo_time)) *
-                        (bb_min[0] + bb_max[0]);
-            center(1) = (0.25 + 0.02 * sin(2 * M_PI * pseudo_time)) *
-                        (bb_min[1] + bb_max[1]);
-            for (int i = 2; i < x.Size(); i++)
-                center(i) = 0.25 * (bb_min[i] + bb_max[i]);
-            rho += amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,
-                                        2.0));
+            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
+            for (int i = 1; i < x.Size(); i++)
+                center(i) = 0.0;
+            rho += amplitude * scaling_factor * std::exp(-50 / scaling_factor *
+                    x.DistanceSquaredTo(center));
         }
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
-            center(0) = (0.75 - 0.02 * cos(2 * M_PI * pseudo_time)) *
-                        (bb_min[0] + bb_max[0]);
-            center(1) = (0.75 - 0.02 * sin(2 * M_PI * pseudo_time)) *
-                        (bb_min[1] + bb_max[1]);
-            for (int i = 2; i < x.Size(); i++)
-                center(i) = 0.75 * (bb_min[i] + bb_max[i]);
-            rho += amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,
-                                        2.0));
+            center(0) = (23 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
+            for (int i = 1; i < x.Size(); i++)
+                center(i) = 0.0;
+            rho += amplitude * scaling_factor * std::exp(-50 / scaling_factor *
+                    x.DistanceSquaredTo(center));
         }
         return rho;
     }

From f6005752f981d1ec911fc41b69701be3e8c44a51 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Mon, 19 Aug 2024 07:28:08 -0700
Subject: [PATCH 36/50] Add variables

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 38 +++++++++++--------
 1 file changed, 23 insertions(+), 15 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 95c0f5562..ee93f91f9 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -70,11 +70,13 @@ using namespace mfem;
 double Conductivity(const Vector &x);
 double Potential(const Vector &x);
 int problem = 1;
-double amplitude = -800.0;
-double pseudo_time = 0.0;
-int potential_well_switch = 0;
-Vector bb_min, bb_max;
-double h_min, h_max, k_min, k_max;
+
+double amplitude = -800.0; // amplitude of coefficient fields
+double pseudo_time = 0.0; // potential parametetr
+Vector bb_min, bb_max; // bounding box
+double h_min, h_max, k_min, k_max; // mesh size
+double alpha = 1.0; // scaling factor
+int potential_well_switch = 0; // switch for separate training
 
 int main(int argc, char *argv[])
 {
@@ -842,8 +844,10 @@ double Conductivity(const Vector &x)
                 return 1.0;
         return 1.0 + amplitude;
     case 3:
-    case 4:
         return 1.0;
+    case 4:
+        alpha = (bb_max[0] - bb_min[0]) / 18.0;
+        return pow(alpha, -2.0);
     }
     return 0.0;
 }
@@ -852,7 +856,8 @@ double Potential(const Vector &x)
 {
     Vector center(x.Size());
     double rho = 0.0;
-    double scaling_factor = pow((bb_max[0] - bb_min[0]) / 18.0, 2.0);
+    double r_sq = 0.0;
+    double radius_sq; // scaling factor
     switch (problem)
     {
     case 1:
@@ -863,26 +868,29 @@ double Potential(const Vector &x)
         center(0) = 0.5 * (bb_min[0] + bb_max[0]) + pseudo_time * h_max;
         for (int i = 1; i < x.Size(); i++)
             center(i) = 0.5 * (bb_min[i] + bb_max[i]);
-        return amplitude * std::exp(-x.DistanceSquaredTo(center) / pow(4.0 * h_max,
-                                    2.0));
+        r_sq = x.DistanceSquaredTo(center);
+        radius_sq = pow(4.0 * h_max, 2.0);
+        return amplitude * std::exp(-r_sq / (2 * radius_sq));
     case 4:
+        alpha = (bb_max[0] - bb_min[0]) / 18.0;
+        radius_sq = pow(0.1 * alpha, 2.0);
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
             center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
             for (int i = 1; i < x.Size(); i++)
-                center(i) = 0.0;
-            rho += amplitude * scaling_factor * std::exp(-50 / scaling_factor *
-                    x.DistanceSquaredTo(center));
+                center(i) = 0.5;
+            r_sq = x.DistanceSquaredTo(center);
+            rho += amplitude * std::exp(-r_sq / (2 * radius_sq));
         }
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
             center(0) = (23 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
             for (int i = 1; i < x.Size(); i++)
-                center(i) = 0.0;
-            rho += amplitude * scaling_factor * std::exp(-50 / scaling_factor *
-                    x.DistanceSquaredTo(center));
+                center(i) = 0.5;
+            r_sq = x.DistanceSquaredTo(center);
+            rho += amplitude * std::exp(-r_sq / (2 * radius_sq));
         }
         return rho;
     }

From 4384d13e2bf9548d5d9cfeceae338290347279d4 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Mon, 19 Aug 2024 08:28:52 -0700
Subject: [PATCH 37/50] Add GTH potential

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 82 +++++++++++++++----
 1 file changed, 67 insertions(+), 15 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index ee93f91f9..c9443d5b6 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -75,7 +75,7 @@ double amplitude = -800.0; // amplitude of coefficient fields
 double pseudo_time = 0.0; // potential parametetr
 Vector bb_min, bb_max; // bounding box
 double h_min, h_max, k_min, k_max; // mesh size
-double alpha = 1.0; // scaling factor
+double L = 1.0; // scaling factor
 int potential_well_switch = 0; // switch for separate training
 
 int main(int argc, char *argv[])
@@ -844,20 +844,21 @@ double Conductivity(const Vector &x)
                 return 1.0;
         return 1.0 + amplitude;
     case 3:
-        return 1.0;
     case 4:
-        alpha = (bb_max[0] - bb_min[0]) / 18.0;
-        return pow(alpha, -2.0);
+        return 1.0;
+    case 5:
+        L = (bb_max[0] - bb_min[0]) / 18.0;
+        return pow(L, -2.0);
     }
-    return 0.0;
+    return 1.0;
 }
 
 double Potential(const Vector &x)
 {
     Vector center(x.Size());
     double rho = 0.0;
-    double r_sq = 0.0;
-    double radius_sq; // scaling factor
+    double d_sq, radius_sq;
+
     switch (problem)
     {
     case 1:
@@ -868,20 +869,21 @@ double Potential(const Vector &x)
         center(0) = 0.5 * (bb_min[0] + bb_max[0]) + pseudo_time * h_max;
         for (int i = 1; i < x.Size(); i++)
             center(i) = 0.5 * (bb_min[i] + bb_max[i]);
-        r_sq = x.DistanceSquaredTo(center);
+        d_sq = x.DistanceSquaredTo(center);
         radius_sq = pow(4.0 * h_max, 2.0);
-        return amplitude * std::exp(-r_sq / (2 * radius_sq));
+        return amplitude * std::exp(-d_sq / (2 * radius_sq));
     case 4:
-        alpha = (bb_max[0] - bb_min[0]) / 18.0;
-        radius_sq = pow(0.1 * alpha, 2.0);
+        L = (bb_max[0] - bb_min[0]) / 18.0;
+        radius_sq = pow(0.1 * L, 2.0);
+
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
             center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
             for (int i = 1; i < x.Size(); i++)
                 center(i) = 0.5;
-            r_sq = x.DistanceSquaredTo(center);
-            rho += amplitude * std::exp(-r_sq / (2 * radius_sq));
+            d_sq = x.DistanceSquaredTo(center);
+            rho += amplitude * std::exp(-d_sq / (2 * radius_sq));
         }
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
@@ -889,8 +891,58 @@ double Potential(const Vector &x)
             center(0) = (23 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
             for (int i = 1; i < x.Size(); i++)
                 center(i) = 0.5;
-            r_sq = x.DistanceSquaredTo(center);
-            rho += amplitude * std::exp(-r_sq / (2 * radius_sq));
+            d_sq = x.DistanceSquaredTo(center);
+            rho += amplitude * std::exp(-d_sq / (2 * radius_sq));
+        }
+        return rho;
+    case 5:
+        L = (bb_max[0] - bb_min[0]) / 18.0;
+        // coefficients for H
+        double rloc = 0.4;
+        double c1 = -14.081155;
+        double c2 = 9.626220;
+        double c3 = -1.783616;
+        double c4 = 0.085152;
+        double zion = 3.0;
+
+        radius_sq = pow(0.1 * L, 2.0);
+        double alpha;
+
+        if (potential_well_switch == 0 || potential_well_switch == 1)
+        {
+            // add well with first center
+            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
+            for (int i = 1; i < x.Size(); i++)
+                center(i) = 0.5;
+            d_sq = x.DistanceSquaredTo(center);
+            alpha = d_sq / pow(L * rloc, 2.0);
+            rho += exp(-0.5 * alpha) * (c1 + c2 * alpha + c3 * alpha * alpha + c4 * alpha * alpha * alpha);
+            if (d_sq > 1.e-16)
+            {
+                rho -= zion * erf(sqrt(d_sq) / (sqrt(2.0) * rloc)) / sqrt(d_sq);
+            }
+            else
+            {
+                rho -= zion * sqrt(2.0) / (sqrt(M_PI) * rloc);
+            }
+        }
+        if (potential_well_switch == 0 || potential_well_switch == 2)
+        {
+            // add well with second center
+            center(0) = (23 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
+            for (int i = 1; i < x.Size(); i++)
+                center(i) = 0.5;
+            d_sq = x.DistanceSquaredTo(center);
+            alpha = d_sq / pow(L * rloc, 2.0);
+            rho += exp(-0.5 * alpha) * (c1 + c2 * alpha + c3 * alpha * alpha + c4 * alpha * alpha * alpha);
+            if (d_sq > 1.e-16)
+            {
+                rho -= zion * erf(sqrt(d_sq) / (sqrt(2.0) * rloc)) / sqrt(d_sq);
+            }
+            else
+            {
+                rho -= zion * sqrt(2.0) / (sqrt(M_PI) * rloc);
+            }
         }
         return rho;
     }

From 7d494bd9a5f5ef58573cd657cf5c1e77dfe789af Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 21 Aug 2024 20:43:14 -0700
Subject: [PATCH 38/50] Slight modification

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 61 +++++++++++++------
 1 file changed, 42 insertions(+), 19 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index c9443d5b6..33770e70e 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -847,6 +847,7 @@ double Conductivity(const Vector &x)
     case 4:
         return 1.0;
     case 5:
+    case 6:
         L = (bb_max[0] - bb_min[0]) / 18.0;
         return pow(L, -2.0);
     }
@@ -865,38 +866,64 @@ double Potential(const Vector &x)
     case 2:
         return 0.0;
     case 3:
+        radius_sq = pow(4.0 * h_max, 2.0);
         // center = (0.5 + t*h, 0.5)
         center(0) = 0.5 * (bb_min[0] + bb_max[0]) + pseudo_time * h_max;
         for (int i = 1; i < x.Size(); i++)
             center(i) = 0.5 * (bb_min[i] + bb_max[i]);
         d_sq = x.DistanceSquaredTo(center);
-        radius_sq = pow(4.0 * h_max, 2.0);
         return amplitude * std::exp(-d_sq / (2 * radius_sq));
     case 4:
-        L = (bb_max[0] - bb_min[0]) / 18.0;
-        radius_sq = pow(0.1 * L, 2.0);
-
+        radius_sq = pow(4.0 * h_max, 2.0);
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
-            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
-            for (int i = 1; i < x.Size(); i++)
-                center(i) = 0.5;
+            center(0) = (0.25 + 0.02 * cos(2 * M_PI * pseudo_time)) *
+                        (bb_min[0] + bb_max[0]);
+            center(1) = (0.25 + 0.02 * sin(2 * M_PI * pseudo_time)) *
+                        (bb_min[1] + bb_max[1]);
+            for (int i = 2; i < x.Size(); i++)
+                center(i) = 0.25 * (bb_min[i] + bb_max[i]);
             d_sq = x.DistanceSquaredTo(center);
             rho += amplitude * std::exp(-d_sq / (2 * radius_sq));
         }
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
-            center(0) = (23 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
-            for (int i = 1; i < x.Size(); i++)
-                center(i) = 0.5;
+            center(0) = (0.75 - 0.02 * cos(2 * M_PI * pseudo_time)) *
+                        (bb_min[0] + bb_max[0]);
+            center(1) = (0.75 - 0.02 * sin(2 * M_PI * pseudo_time)) *
+                        (bb_min[1] + bb_max[1]);
+            for (int i = 2; i < x.Size(); i++)
+                center(i) = 0.25 * (bb_min[i] + bb_max[i]);
             d_sq = x.DistanceSquaredTo(center);
             rho += amplitude * std::exp(-d_sq / (2 * radius_sq));
         }
         return rho;
     case 5:
+        radius_sq = pow(0.1, 2.0);
+        for (int i = 1; i < x.Size(); i++)
+            center(i) = (bb_min[i] + bb_max[i]) / 2;
+        if (potential_well_switch == 0 || potential_well_switch == 1)
+        {
+            // add well with first center
+            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
+            d_sq = x.DistanceSquaredTo(center);
+            rho += amplitude * std::exp(-d_sq / (2 * radius_sq));
+        }
+        if (potential_well_switch == 0 || potential_well_switch == 2)
+        {
+            // add well with second center
+            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
+            d_sq = x.DistanceSquaredTo(center);
+            rho += amplitude * std::exp(-d_sq / (2 * radius_sq));
+        }
+        return rho;
+    case 6:
         L = (bb_max[0] - bb_min[0]) / 18.0;
+        for (int i = 1; i < x.Size(); i++)
+            center(i) = (bb_min[i] + bb_max[i]) / 2;
+
         // coefficients for H
         double rloc = 0.4;
         double c1 = -14.081155;
@@ -905,18 +932,15 @@ double Potential(const Vector &x)
         double c4 = 0.085152;
         double zion = 3.0;
 
-        radius_sq = pow(0.1 * L, 2.0);
         double alpha;
-
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
             center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
-            for (int i = 1; i < x.Size(); i++)
-                center(i) = 0.5;
             d_sq = x.DistanceSquaredTo(center);
             alpha = d_sq / pow(L * rloc, 2.0);
-            rho += exp(-0.5 * alpha) * (c1 + c2 * alpha + c3 * alpha * alpha + c4 * alpha * alpha * alpha);
+            rho += exp(-0.5 * alpha) * (c1 + c2 * alpha + c3 * alpha * alpha + c4 * alpha *
+                                        alpha * alpha);
             if (d_sq > 1.e-16)
             {
                 rho -= zion * erf(sqrt(d_sq) / (sqrt(2.0) * rloc)) / sqrt(d_sq);
@@ -929,12 +953,11 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
-            center(0) = (23 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
-            for (int i = 1; i < x.Size(); i++)
-                center(i) = 0.5;
+            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
             d_sq = x.DistanceSquaredTo(center);
             alpha = d_sq / pow(L * rloc, 2.0);
-            rho += exp(-0.5 * alpha) * (c1 + c2 * alpha + c3 * alpha * alpha + c4 * alpha * alpha * alpha);
+            rho += exp(-0.5 * alpha) * (c1 + c2 * alpha + c3 * alpha * alpha + c4 * alpha *
+                                        alpha * alpha);
             if (d_sq > 1.e-16)
             {
                 rho -= zion * erf(sqrt(d_sq) / (sqrt(2.0) * rloc)) / sqrt(d_sq);

From f466b6f5b11bd2da3039cb12cd563cbf143218cd Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 21 Aug 2024 21:13:43 -0700
Subject: [PATCH 39/50] Slight updates

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 33770e70e..22c322f92 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -849,7 +849,7 @@ double Conductivity(const Vector &x)
     case 5:
     case 6:
         L = (bb_max[0] - bb_min[0]) / 18.0;
-        return pow(L, -2.0);
+        return pow(L, 2.0);
     }
     return 1.0;
 }
@@ -895,15 +895,16 @@ double Potential(const Vector &x)
             center(1) = (0.75 - 0.02 * sin(2 * M_PI * pseudo_time)) *
                         (bb_min[1] + bb_max[1]);
             for (int i = 2; i < x.Size(); i++)
-                center(i) = 0.25 * (bb_min[i] + bb_max[i]);
+                center(i) = 0.75 * (bb_min[i] + bb_max[i]);
             d_sq = x.DistanceSquaredTo(center);
             rho += amplitude * std::exp(-d_sq / (2 * radius_sq));
         }
         return rho;
     case 5:
-        radius_sq = pow(0.1, 2.0);
+        L = (bb_max[0] - bb_min[0]) / 18.0;
         for (int i = 1; i < x.Size(); i++)
             center(i) = (bb_min[i] + bb_max[i]) / 2;
+        radius_sq = pow(0.2 * L, 2.0);
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center

From 7441b0f8cde09d5234e465c51f1d14763af7d112 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 22 Aug 2024 12:54:38 -0700
Subject: [PATCH 40/50] superposing gaussians

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 22c322f92..ac2ce7f2f 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -858,6 +858,7 @@ double Potential(const Vector &x)
 {
     Vector center(x.Size());
     double rho = 0.0;
+    double rho_bg = -5.0;
     double d_sq, radius_sq;
 
     switch (problem)
@@ -904,20 +905,25 @@ double Potential(const Vector &x)
         L = (bb_max[0] - bb_min[0]) / 18.0;
         for (int i = 1; i < x.Size(); i++)
             center(i) = (bb_min[i] + bb_max[i]) / 2;
-        radius_sq = pow(0.2 * L, 2.0);
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
             center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
             d_sq = x.DistanceSquaredTo(center);
-            rho += amplitude * std::exp(-d_sq / (2 * radius_sq));
+            radius_sq = pow(0.28 * L, 2.0);
+            rho += -29.0 * std::exp(-d_sq / (2 * radius_sq));
+            radius_sq = pow(1.08 * L, 2.0);
+            rho += -3.0 * std::exp(-d_sq / (2 * radius_sq));
         }
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
             center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
             d_sq = x.DistanceSquaredTo(center);
-            rho += amplitude * std::exp(-d_sq / (2 * radius_sq));
+            radius_sq = pow(0.28 * L, 2.0);
+            rho += -29.0 * std::exp(-d_sq / (2 * radius_sq));
+            radius_sq = pow(1.08 * L, 2.0);
+            rho += -3.0 * std::exp(-d_sq / (2 * radius_sq));
         }
         return rho;
     case 6:

From 77516d43f2f2777349b4f8ca525bcb59cbbb84cc Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Wed, 28 Aug 2024 20:57:53 -0700
Subject: [PATCH 41/50] Minor changes to parameters

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index ac2ce7f2f..00f109ae1 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -911,9 +911,9 @@ double Potential(const Vector &x)
             center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
             d_sq = x.DistanceSquaredTo(center);
             radius_sq = pow(0.28 * L, 2.0);
-            rho += -29.0 * std::exp(-d_sq / (2 * radius_sq));
+            rho += -28.9 * std::exp(-d_sq / (2 * radius_sq));
             radius_sq = pow(1.08 * L, 2.0);
-            rho += -3.0 * std::exp(-d_sq / (2 * radius_sq));
+            rho += -3.6 * std::exp(-d_sq / (2 * radius_sq));
         }
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
@@ -921,9 +921,9 @@ double Potential(const Vector &x)
             center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
             d_sq = x.DistanceSquaredTo(center);
             radius_sq = pow(0.28 * L, 2.0);
-            rho += -29.0 * std::exp(-d_sq / (2 * radius_sq));
+            rho += -28.9 * std::exp(-d_sq / (2 * radius_sq));
             radius_sq = pow(1.08 * L, 2.0);
-            rho += -3.0 * std::exp(-d_sq / (2 * radius_sq));
+            rho += -3.6 * std::exp(-d_sq / (2 * radius_sq));
         }
         return rho;
     case 6:

From 67264360ef0c51f0e27624d03d563fa63cc9924d Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 29 Aug 2024 09:27:31 -0700
Subject: [PATCH 42/50] Minor changes

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 00f109ae1..0e2f65a96 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -908,7 +908,8 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
-            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
+            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 -
+                        0.5 * h_max * pseudo_time * (bb_min[0] + bb_max[0]);
             d_sq = x.DistanceSquaredTo(center);
             radius_sq = pow(0.28 * L, 2.0);
             rho += -28.9 * std::exp(-d_sq / (2 * radius_sq));
@@ -918,7 +919,8 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
-            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
+            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 +
+                        0.5 * h_max * pseudo_time * (bb_min[0] + bb_max[0]);
             d_sq = x.DistanceSquaredTo(center);
             radius_sq = pow(0.28 * L, 2.0);
             rho += -28.9 * std::exp(-d_sq / (2 * radius_sq));

From 29631ce0e47b6ef98c914dad735bff8299ceb40a Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 5 Sep 2024 08:14:51 -0700
Subject: [PATCH 43/50] Add loading

---
 .../prom/elliptic_eigenproblem_global_rom.cpp  | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 0e2f65a96..36eff427b 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -458,6 +458,24 @@ int main(int argc, char *argv[])
         lobpcg->SetMassMatrix(*M);
         lobpcg->SetOperator(*A);
 
+        if (id > 0)
+        {
+            std::string snapshot_filename = baseName + "par0_snapshot";
+            CAROM::BasisGenerator* aux_generator = new CAROM::BasisGenerator(*options, isIncremental, basis_filename);
+            aux_generator->loadSamples(snapshot_filename, "snapshot");
+            const CAROM::Matrix* snapshot_mat_carom = aux_generator->getSnapshotMatrix();
+            CAROM::Vector* snapshot_vec_carom = nullptr;
+            Vector** snapshot_vecs = new Vector*[nev];
+            for (int j = 0; j < nev; j++)
+            {
+                snapshot_vecs[j] = new Vector(size);
+                snapshot_mat_carom->getColumn(j, snapshot_vec_carom);
+                snapshot_vecs[j]->SetData(snapshot_vec_carom->getData());
+            }
+            delete aux_generator;
+            lobpcg->SetInitialVectors(nev, snapshot_vecs);
+        }
+
         // 13. Compute the eigenmodes and extract the array of eigenvalues. Define a
         //     parallel grid function to represent each of the eigenmodes returned by
         //     the solver.

From 6de0d5a31c9f08785e25e54d73bc7faab760b9b1 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 5 Sep 2024 08:49:38 -0700
Subject: [PATCH 44/50] Use HypreParVector I/O instead of libROM I/O for this
 task

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 23 ++++++++++---------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 36eff427b..5c03ed2c3 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -460,19 +460,13 @@ int main(int argc, char *argv[])
 
         if (id > 0)
         {
-            std::string snapshot_filename = baseName + "par0_snapshot";
-            CAROM::BasisGenerator* aux_generator = new CAROM::BasisGenerator(*options, isIncremental, basis_filename);
-            aux_generator->loadSamples(snapshot_filename, "snapshot");
-            const CAROM::Matrix* snapshot_mat_carom = aux_generator->getSnapshotMatrix();
-            CAROM::Vector* snapshot_vec_carom = nullptr;
-            Vector** snapshot_vecs = new Vector*[nev];
-            for (int j = 0; j < nev; j++)
+            HypreParVector** snapshot_vecs = new HypreParVector*[nev];
+            for (int i = 0; i < nev; i++)
             {
-                snapshot_vecs[j] = new Vector(size);
-                snapshot_mat_carom->getColumn(j, snapshot_vec_carom);
-                snapshot_vecs[j]->SetData(snapshot_vec_carom->getData());
+                std::string snapshot_filename = baseName + "ref_snapshot_" + std::to_string(i);
+                snapshot_vecs[i] = new HypreParVector();
+                snapshot_vecs[i]->Read(MPI_COMM_WORLD, snapshot_filename.c_str());
             }
-            delete aux_generator;
             lobpcg->SetInitialVectors(nev, snapshot_vecs);
         }
 
@@ -487,6 +481,13 @@ int main(int argc, char *argv[])
         // 14. take and write snapshots for ROM
         for (int i = 0; i < nev; i++)
         {
+            if (id == 0)
+            {
+                std::string snapshot_filename = baseName + "ref_snapshot_" + std::to_string(i);
+               const HypreParVector snapshot_vec = lobpcg->GetEigenvector(i);
+               snapshot_vec.Print(snapshot_filename.c_str());
+            }
+
             if (myid == 0)
             {
                 std::cout << " Eigenvalue " << i << ": " << eigenvalues[i] << "\n";

From 407227b1e439551696b6e3571ad63f3b9a0ab554 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 5 Sep 2024 09:06:25 -0700
Subject: [PATCH 45/50] Astyle

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 5c03ed2c3..3ae8dddd3 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -484,8 +484,8 @@ int main(int argc, char *argv[])
             if (id == 0)
             {
                 std::string snapshot_filename = baseName + "ref_snapshot_" + std::to_string(i);
-               const HypreParVector snapshot_vec = lobpcg->GetEigenvector(i);
-               snapshot_vec.Print(snapshot_filename.c_str());
+                const HypreParVector snapshot_vec = lobpcg->GetEigenvector(i);
+                snapshot_vec.Print(snapshot_filename.c_str());
             }
 
             if (myid == 0)

From 7adae585e072214e0ef255d78f00340e5cdc0696 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 12 Sep 2024 08:31:38 -0700
Subject: [PATCH 46/50] Add option to compare initialization

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 58 ++++++++++++-------
 1 file changed, 38 insertions(+), 20 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 3ae8dddd3..05df99a62 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -111,6 +111,7 @@ int main(int argc, char *argv[])
     double ef = 1.0;
     int rdim = -1;
     int verbose_level = 0;
+    bool prescribe_init = false;
 
     int precision = 8;
     cout.precision(precision);
@@ -163,6 +164,8 @@ int main(int argc, char *argv[])
                    "Reduced dimension.");
     args.AddOption(&verbose_level, "-v", "--verbose",
                    "Set the verbosity level of the LOBPCG solver and preconditioner. 0 is off.");
+    args.AddOption(&prescribe_init, "-init", "--init", "-no-init", "--no-init",
+                   "Prescribe custom Initialization of LOBPCG.");
     args.AddOption(&lobpcg_niter, "-fi", "--fom-iter",
                    "Number of iterations for the LOBPCG solver.");
 #ifdef MFEM_USE_SUPERLU
@@ -458,16 +461,30 @@ int main(int argc, char *argv[])
         lobpcg->SetMassMatrix(*M);
         lobpcg->SetOperator(*A);
 
-        if (id > 0)
+        if (prescribe_init && (fom || (offline && id > 0)))
         {
             HypreParVector** snapshot_vecs = new HypreParVector*[nev];
             for (int i = 0; i < nev; i++)
             {
                 std::string snapshot_filename = baseName + "ref_snapshot_" + std::to_string(i);
-                snapshot_vecs[i] = new HypreParVector();
-                snapshot_vecs[i]->Read(MPI_COMM_WORLD, snapshot_filename.c_str());
+                std::ifstream snapshot_infile(snapshot_filename + ".0");
+                if (!snapshot_infile.good())
+                {
+                    if (myid == 0) std::cout << "Failed to load " << snapshot_filename << std::endl;
+                    break;
+                }
+                else
+                {
+                    snapshot_vecs[i] = new HypreParVector();
+                    snapshot_vecs[i]->Read(MPI_COMM_WORLD, snapshot_filename.c_str());
+                    if (myid == 0) std::cout << "Loaded " << snapshot_filename << std::endl;
+                    if (i == nev - 1)
+                    {
+                        lobpcg->SetInitialVectors(nev, snapshot_vecs);
+                        if (myid == 0) std::cout << "LOBPCG initial vectors set" << std::endl;
+                    }
+                }
             }
-            lobpcg->SetInitialVectors(nev, snapshot_vecs);
         }
 
         // 13. Compute the eigenmodes and extract the array of eigenvalues. Define a
@@ -481,13 +498,6 @@ int main(int argc, char *argv[])
         // 14. take and write snapshots for ROM
         for (int i = 0; i < nev; i++)
         {
-            if (id == 0)
-            {
-                std::string snapshot_filename = baseName + "ref_snapshot_" + std::to_string(i);
-                const HypreParVector snapshot_vec = lobpcg->GetEigenvector(i);
-                snapshot_vec.Print(snapshot_filename.c_str());
-            }
-
             if (myid == 0)
             {
                 std::cout << " Eigenvalue " << i << ": " << eigenvalues[i] << "\n";
@@ -497,6 +507,12 @@ int main(int argc, char *argv[])
                 eigenfunction_i = lobpcg->GetEigenvector(i);
                 eigenfunction_i /= sqrt(InnerProduct(eigenfunction_i, eigenfunction_i));
                 generator->takeSample(eigenfunction_i.GetData());
+                if (prescribe_init && id == 0)
+                {
+                    std::string snapshot_filename = baseName + "ref_snapshot_" + std::to_string(i);
+                    const HypreParVector snapshot_vec = lobpcg->GetEigenvector(i);
+                    snapshot_vec.Print(snapshot_filename.c_str());
+                }
             }
         }
 
@@ -604,6 +620,9 @@ int main(int argc, char *argv[])
         // Calculate errors of eigenvectors
         ostringstream mode_name_fom;
         Vector mode_fom(numRowRB);
+        double eigenvalue_threshold = max(1e-6,
+                                          0.02 * (eigenvalues_fom[nev-1] - eigenvalues_fom[0]));
+        std::cout << "eigenvalue_threshold = " << eigenvalue_threshold << std::endl;
         for (int i = 0; i < nev; i++)
         {
             mode_name_fom << "eigenfunction_fom_" << setfill('0') << setw(2) << i << "."
@@ -617,12 +636,12 @@ int main(int argc, char *argv[])
             eigenvector_i = 0.0;
             for (int j = 0; j < nev; j++)
             {
-                if (myid == 0 && verbose_level > 0)
+                if (myid == 0)
                 {
                     std::cout << "correlation_matrix(" << j+1 << "," << i+1 << ") = " <<
                               InnerProduct(mode_fom, modes_rom[j]) << ";" << std::endl;
                 }
-                if (abs(eigenvalues_fom[j] - eigenvalues_fom[i]) < 1e-6)
+                if (abs(eigenvalues_fom[j] - eigenvalues_fom[i]) < eigenvalue_threshold)
                 {
                     eigenvector_i.Add(InnerProduct(mode_fom, modes_rom[j]), modes_rom[j]);
                 }
@@ -927,8 +946,7 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
-            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 -
-                        0.5 * h_max * pseudo_time * (bb_min[0] + bb_max[0]);
+            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - h_max * pseudo_time;
             d_sq = x.DistanceSquaredTo(center);
             radius_sq = pow(0.28 * L, 2.0);
             rho += -28.9 * std::exp(-d_sq / (2 * radius_sq));
@@ -938,8 +956,7 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
-            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 +
-                        0.5 * h_max * pseudo_time * (bb_min[0] + bb_max[0]);
+            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 + h_max * pseudo_time;
             d_sq = x.DistanceSquaredTo(center);
             radius_sq = pow(0.28 * L, 2.0);
             rho += -28.9 * std::exp(-d_sq / (2 * radius_sq));
@@ -951,7 +968,6 @@ double Potential(const Vector &x)
         L = (bb_max[0] - bb_min[0]) / 18.0;
         for (int i = 1; i < x.Size(); i++)
             center(i) = (bb_min[i] + bb_max[i]) / 2;
-
         // coefficients for H
         double rloc = 0.4;
         double c1 = -14.081155;
@@ -964,7 +980,8 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
-            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - pseudo_time * h_max;
+            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 -
+                        h_max * pseudo_time * (bb_min[0] + bb_max[0]) / 2;
             d_sq = x.DistanceSquaredTo(center);
             alpha = d_sq / pow(L * rloc, 2.0);
             rho += exp(-0.5 * alpha) * (c1 + c2 * alpha + c3 * alpha * alpha + c4 * alpha *
@@ -981,7 +998,8 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
-            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 + pseudo_time * h_max;
+            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 +
+                        h_max * pseudo_time * (bb_min[0] + bb_max[0]) / 2;
             d_sq = x.DistanceSquaredTo(center);
             alpha = d_sq / pow(L * rloc, 2.0);
             rho += exp(-0.5 * alpha) * (c1 + c2 * alpha + c3 * alpha * alpha + c4 * alpha *

From f75855edb4208f10fc925a24084e30e3ee536457 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 12 Sep 2024 09:02:22 -0700
Subject: [PATCH 47/50] Minor change to example

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 6a86aa186..0ba5a72d0 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -946,8 +946,7 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
-            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 -
-                        0.5 * h_max * pseudo_time * (bb_min[0] + bb_max[0]);
+            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - h_max * pseudo_time;
             d_sq = x.DistanceSquaredTo(center);
             radius_sq = pow(0.28 * L, 2.0);
             rho += -28.9 * std::exp(-d_sq / (2 * radius_sq));
@@ -957,8 +956,7 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
-            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 +
-                        0.5 * h_max * pseudo_time * (bb_min[0] + bb_max[0]);
+            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 + h_max * pseudo_time;
             d_sq = x.DistanceSquaredTo(center);
             radius_sq = pow(0.28 * L, 2.0);
             rho += -28.9 * std::exp(-d_sq / (2 * radius_sq));
@@ -983,8 +981,7 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 1)
         {
             // add well with first center
-            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 -
-                        h_max * pseudo_time * (bb_min[0] + bb_max[0]) / 2;
+            center(0) = (23 * bb_min[0] + 13 * bb_max[0]) / 36 - h_max * pseudo_time;
             d_sq = x.DistanceSquaredTo(center);
             alpha = d_sq / pow(L * rloc, 2.0);
             rho += exp(-0.5 * alpha) * (c1 + c2 * alpha + c3 * alpha * alpha + c4 * alpha *
@@ -1001,8 +998,7 @@ double Potential(const Vector &x)
         if (potential_well_switch == 0 || potential_well_switch == 2)
         {
             // add well with second center
-            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 +
-                        h_max * pseudo_time * (bb_min[0] + bb_max[0]) / 2;
+            center(0) = (13 * bb_min[0] + 23 * bb_max[0]) / 36 + h_max * pseudo_time;
             d_sq = x.DistanceSquaredTo(center);
             alpha = d_sq / pow(L * rloc, 2.0);
             rho += exp(-0.5 * alpha) * (c1 + c2 * alpha + c3 * alpha * alpha + c4 * alpha *

From cd9dda8c84463c10a17131de35a6991527f1da12 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 12 Sep 2024 11:42:53 -0700
Subject: [PATCH 48/50] Abort if no files found

---
 .../prom/elliptic_eigenproblem_global_rom.cpp | 24 +++++++------------
 1 file changed, 8 insertions(+), 16 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 0ba5a72d0..ab355b991 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -467,24 +467,15 @@ int main(int argc, char *argv[])
             for (int i = 0; i < nev; i++)
             {
                 std::string snapshot_filename = baseName + "ref_snapshot_" + std::to_string(i);
+                if (myid == 0) std::cout << "Loading " << snapshot_filename << std::endl;
                 std::ifstream snapshot_infile(snapshot_filename + ".0");
-                if (!snapshot_infile.good())
-                {
-                    if (myid == 0) std::cout << "Failed to load " << snapshot_filename << std::endl;
-                    break;
-                }
-                else
-                {
-                    snapshot_vecs[i] = new HypreParVector();
-                    snapshot_vecs[i]->Read(MPI_COMM_WORLD, snapshot_filename.c_str());
-                    if (myid == 0) std::cout << "Loaded " << snapshot_filename << std::endl;
-                    if (i == nev - 1)
-                    {
-                        lobpcg->SetInitialVectors(nev, snapshot_vecs);
-                        if (myid == 0) std::cout << "LOBPCG initial vectors set" << std::endl;
-                    }
-                }
+                MFEM_VERIFY(snapshot_infile.good(), "Failed to load vector for initialization");
+                snapshot_vecs[i] = new HypreParVector();
+                snapshot_vecs[i]->Read(MPI_COMM_WORLD, snapshot_filename.c_str());
+                if (myid == 0) std::cout << "Loaded " << snapshot_filename << std::endl;
             }
+            lobpcg->SetInitialVectors(nev, snapshot_vecs);
+            if (myid == 0) std::cout << "LOBPCG initial vectors set" << std::endl;
         }
 
         // 13. Compute the eigenmodes and extract the array of eigenvalues. Define a
@@ -512,6 +503,7 @@ int main(int argc, char *argv[])
                     std::string snapshot_filename = baseName + "ref_snapshot_" + std::to_string(i);
                     const HypreParVector snapshot_vec = lobpcg->GetEigenvector(i);
                     snapshot_vec.Print(snapshot_filename.c_str());
+                    if (myid == 0) std::cout << "Saved " << snapshot_filename << std::endl;
                 }
             }
         }

From 1222a809090512e4941d6d42fd8084a2685a053d Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 19 Sep 2024 07:03:02 -0700
Subject: [PATCH 49/50] Minor change

---
 examples/prom/elliptic_eigenproblem_global_rom.cpp | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index ab355b991..7c231bd66 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -454,8 +454,8 @@ int main(int argc, char *argv[])
         lobpcg->SetNumModes(nev);
         lobpcg->SetRandomSeed(seed);
         lobpcg->SetPreconditioner(*precond);
-        lobpcg->SetMaxIter(200);
-        lobpcg->SetTol(1e-8);
+        lobpcg->SetMaxIter(lobpcg_niter);
+        lobpcg->SetTol(1e-6);
         lobpcg->SetPrecondUsageMode(1);
         lobpcg->SetPrintLevel(verbose_level);
         lobpcg->SetMassMatrix(*M);
@@ -467,8 +467,7 @@ int main(int argc, char *argv[])
             for (int i = 0; i < nev; i++)
             {
                 std::string snapshot_filename = baseName + "ref_snapshot_" + std::to_string(i);
-                if (myid == 0) std::cout << "Loading " << snapshot_filename << std::endl;
-                std::ifstream snapshot_infile(snapshot_filename + ".0");
+                std::ifstream snapshot_infile(snapshot_filename + "." + std::to_string(myid));
                 MFEM_VERIFY(snapshot_infile.good(), "Failed to load vector for initialization");
                 snapshot_vecs[i] = new HypreParVector();
                 snapshot_vecs[i]->Read(MPI_COMM_WORLD, snapshot_filename.c_str());
@@ -612,8 +611,8 @@ int main(int argc, char *argv[])
         // Calculate errors of eigenvectors
         ostringstream mode_name_fom;
         Vector mode_fom(numRowRB);
-        double eigenvalue_threshold = max(1e-6,
-                                          0.02 * (eigenvalues_fom[nev-1] - eigenvalues_fom[0]));
+        const double eigenvalue_threshold = max(1e-6,
+                                                0.02 * (eigenvalues_fom[nev-1] - eigenvalues_fom[0]));
         std::cout << "eigenvalue_threshold = " << eigenvalue_threshold << std::endl;
         for (int i = 0; i < nev; i++)
         {

From 32d545f52b8f92ee027b9333f0e7145f275395f8 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 26 Sep 2024 08:56:43 -0700
Subject: [PATCH 50/50] Make tolerences user-defined options

---
 .../prom/elliptic_eigenproblem_global_rom.cpp   | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp
index 7c231bd66..815a493bf 100644
--- a/examples/prom/elliptic_eigenproblem_global_rom.cpp
+++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp
@@ -94,10 +94,14 @@ int main(int argc, char *argv[])
     int order = 1;
     int nev = 4;
     int seed = 75;
+    bool prescribe_init = false;
     int lobpcg_niter = 200;
+    double lobpcg_tol = 1e-8;
+    double eig_tol = 1e-6;
     bool slu_solver = false;
     bool sp_solver = false;
     bool cpardiso_solver = false;
+    int verbose_level = 0;
 
     bool visualization = true;
     bool visit = false;
@@ -110,8 +114,6 @@ int main(int argc, char *argv[])
     int nsets = 0;
     double ef = 1.0;
     int rdim = -1;
-    int verbose_level = 0;
-    bool prescribe_init = false;
 
     int precision = 8;
     cout.precision(precision);
@@ -168,6 +170,10 @@ int main(int argc, char *argv[])
                    "Prescribe custom Initialization of LOBPCG.");
     args.AddOption(&lobpcg_niter, "-fi", "--fom-iter",
                    "Number of iterations for the LOBPCG solver.");
+    args.AddOption(&lobpcg_tol, "-tol", "--fom-tol",
+                   "Tolerance for the LOBPCG solver.");
+    args.AddOption(&eig_tol, "-tol", "--fom-tol",
+                   "Tolerance for eigenvalues to be considered equal.");
 #ifdef MFEM_USE_SUPERLU
     args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
                    "--no-superlu", "Use the SuperLU Solver.");
@@ -455,7 +461,7 @@ int main(int argc, char *argv[])
         lobpcg->SetRandomSeed(seed);
         lobpcg->SetPreconditioner(*precond);
         lobpcg->SetMaxIter(lobpcg_niter);
-        lobpcg->SetTol(1e-6);
+        lobpcg->SetTol(lobpcg_tol);
         lobpcg->SetPrecondUsageMode(1);
         lobpcg->SetPrintLevel(verbose_level);
         lobpcg->SetMassMatrix(*M);
@@ -611,9 +617,6 @@ int main(int argc, char *argv[])
         // Calculate errors of eigenvectors
         ostringstream mode_name_fom;
         Vector mode_fom(numRowRB);
-        const double eigenvalue_threshold = max(1e-6,
-                                                0.02 * (eigenvalues_fom[nev-1] - eigenvalues_fom[0]));
-        std::cout << "eigenvalue_threshold = " << eigenvalue_threshold << std::endl;
         for (int i = 0; i < nev; i++)
         {
             mode_name_fom << "eigenfunction_fom_" << setfill('0') << setw(2) << i << "."
@@ -632,7 +635,7 @@ int main(int argc, char *argv[])
                     std::cout << "correlation_matrix(" << j+1 << "," << i+1 << ") = " <<
                               InnerProduct(mode_fom, modes_rom[j]) << ";" << std::endl;
                 }
-                if (abs(eigenvalues_fom[j] - eigenvalues_fom[i]) < eigenvalue_threshold)
+                if (abs(eigenvalues_fom[j] - eigenvalues_fom[i]) < eig_tol)
                 {
                     eigenvector_i.Add(InnerProduct(mode_fom, modes_rom[j]), modes_rom[j]);
                 }