DMD with shifted data

examples/dmd/dg_advection.cpp

@@ -119,7 +119,6 @@ class AIR_prec : public Solver
 };
 #endif
 
-
 class DG_Solver : public Solver
 {
 private:
@@ -195,7 +194,6 @@ class DG_Solver : public Solver
     }
 };
 
-
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
     and advection matrices, and b describes the flow on the boundary. This can
@@ -222,7 +220,6 @@ class FE_Evolution : public TimeDependentOperator
     virtual ~FE_Evolution();
 };
 
-
 int main(int argc, char *argv[])
 {
     // 1. Initialize MPI.
@@ -465,7 +462,6 @@ int main(int argc, char *argv[])
     m->Finalize();
     k->Finalize(skip_zeros);
 
-
     HypreParVector *B = b->ParallelAssemble();
 
     // 9. Define the initial conditions, save the corresponding grid function to
@@ -798,7 +794,6 @@ int main(int argc, char *argv[])
     return 0;
 }
 
-
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
                            const Vector &b_, PrecType prec_type)
@@ -869,7 +864,6 @@ FE_Evolution::~FE_Evolution()
     delete dg_solver;
 }
 
-
 // Velocity coefficient
 void velocity_function(const Vector &x, Vector &v)
 {

examples/dmd/dg_euler.hpp

@@ -56,7 +56,6 @@ class RiemannSolver
                 const Vector &nor, Vector &flux);
 };
 
-
 // Constant (in time) mixed bilinear form multiplying the flux grid function.
 // The form is (vec(v), grad(w)) where the trial space = vector L2 space (mesh
 // dim) and test space = scalar L2 space.

examples/dmd/nonlinear_elasticity.cpp

@@ -176,7 +176,6 @@ class ReducedSystemOperator : public Operator
     virtual ~ReducedSystemOperator();
 };
 
-
 /** Function representing the elastic energy density for the given hyperelastic
     model+deformation. Used in HyperelasticOperator::GetElasticEnergyDensity. */
 class ElasticEnergyCoefficient : public Coefficient
@@ -201,7 +200,6 @@ void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
                ParGridFunction *field, const char *field_name = NULL,
                bool init_vis = false);
 
-
 int main(int argc, char *argv[])
 {
     // 1. Initialize MPI.
@@ -606,11 +604,6 @@ int main(int argc, char *argv[])
         std::cout << "Both rdim and ef are set. ef will be ignored." << std::endl;
     }
 
-    dmd_training_timer.Start();
-
-
-    dmd_training_timer.Stop();
-
     Vector true_solution_x(x_gf.GetTrueVector().Size());
     true_solution_x = x_gf.GetTrueVector();
 
@@ -767,7 +760,6 @@ void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
     out << flush;
 }
 
-
 ReducedSystemOperator::ReducedSystemOperator(
     ParBilinearForm *M_, ParBilinearForm *S_, ParNonlinearForm *H_,
     const Array<int> &ess_tdof_list_)
@@ -814,7 +806,6 @@ ReducedSystemOperator::~ReducedSystemOperator()
     delete Jacobian;
 }
 
-
 HyperelasticOperator::HyperelasticOperator(ParFiniteElementSpace &f,
         Array<int> &ess_bdr, double visc,
         double mu, double K)
@@ -952,7 +943,6 @@ HyperelasticOperator::~HyperelasticOperator()
     delete Mmat;
 }
 
-
 double ElasticEnergyCoefficient::Eval(ElementTransformation &T,
                                       const IntegrationPoint &ip)
 {
@@ -962,7 +952,6 @@ double ElasticEnergyCoefficient::Eval(ElementTransformation &T,
     return model.EvalW(J)/J.Det(); // in deformed configuration
 }
 
-
 void InitialDeformation(const Vector &x, Vector &y)
 {
     // set the initial configuration to be the same as the reference, stress

examples/prom/dg_advection_global_rom.cpp

@@ -146,7 +146,6 @@ class AIR_prec : public Solver
 };
 #endif
 
-
 class DG_Solver : public Solver
 {
 private:
@@ -222,7 +221,6 @@ class DG_Solver : public Solver
     }
 };
 
-
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
     and advection matrices, and b describes the flow on the boundary. This can
@@ -245,7 +243,6 @@ class ROM_FE_Evolution : public TimeDependentOperator
     virtual ~ROM_FE_Evolution();
 };
 
-
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
     and advection matrices, and b describes the flow on the boundary. This can
@@ -272,7 +269,6 @@ class FE_Evolution : public TimeDependentOperator
     virtual ~FE_Evolution();
 };
 
-
 int main(int argc, char *argv[])
 {
     // 1. Initialize MPI.
@@ -548,7 +544,6 @@ int main(int argc, char *argv[])
     m->Finalize();
     k->Finalize(skip_zeros);
 
-
     HypreParVector *B = b->ParallelAssemble();
 
     // 9. Define the initial conditions, save the corresponding grid function to
@@ -955,7 +950,6 @@ int main(int argc, char *argv[])
     return 0;
 }
 
-
 // Implementation of class ROM_FE_Evolution
 ROM_FE_Evolution::ROM_FE_Evolution(DenseMatrix* M_, DenseMatrix* K_, Vector* b_,
                                    Vector* u_init_hat_, int num_cols)
@@ -1011,7 +1005,6 @@ ROM_FE_Evolution::~ROM_FE_Evolution()
     delete A_inv;
 }
 
-
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
                            const Vector &b_, PrecType prec_type)
@@ -1082,7 +1075,6 @@ FE_Evolution::~FE_Evolution()
     delete dg_solver;
 }
 
-
 // Velocity coefficient
 void velocity_function(const Vector &x, Vector &v)
 {

examples/prom/dg_advection_local_rom_matrix_interp.cpp

@@ -158,7 +158,6 @@ class AIR_prec : public Solver
 };
 #endif
 
-
 class DG_Solver : public Solver
 {
 private:
@@ -234,7 +233,6 @@ class DG_Solver : public Solver
     }
 };
 
-
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
     and advection matrices, and b describes the flow on the boundary. This can
@@ -257,7 +255,6 @@ class ROM_FE_Evolution : public TimeDependentOperator
     virtual ~ROM_FE_Evolution();
 };
 
-
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
     and advection matrices, and b describes the flow on the boundary. This can
@@ -284,7 +281,6 @@ class FE_Evolution : public TimeDependentOperator
     virtual ~FE_Evolution();
 };
 
-
 int main(int argc, char *argv[])
 {
     // 1. Initialize MPI.
@@ -570,7 +566,6 @@ int main(int argc, char *argv[])
     m->Finalize();
     k->Finalize(skip_zeros);
 
-
     HypreParVector *B = b->ParallelAssemble();
 
     // 9. Define the initial conditions, save the corresponding grid function to
@@ -1065,7 +1060,6 @@ int main(int argc, char *argv[])
     return 0;
 }
 
-
 // Implementation of class ROM_FE_Evolution
 ROM_FE_Evolution::ROM_FE_Evolution(DenseMatrix* M_, DenseMatrix* K_, Vector* b_,
                                    Vector* u_init_hat_, int num_cols)
@@ -1121,7 +1115,6 @@ ROM_FE_Evolution::~ROM_FE_Evolution()
     delete A_inv;
 }
 
-
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
                            const Vector &b_, PrecType prec_type)
@@ -1192,7 +1185,6 @@ FE_Evolution::~FE_Evolution()
     delete dg_solver;
 }
 
-
 // Velocity coefficient
 void velocity_function(const Vector &x, Vector &v)
 {

examples/prom/mixed_nonlinear_diffusion.cpp

@@ -67,7 +67,6 @@
 #include "hyperreduction/S_OPT.h"
 #include "mfem/SampleMesh.hpp"
 
-
 typedef enum {ANALYTIC, INIT_STEP} PROBLEM;
 
 typedef enum {RSPACE, WSPACE} FESPACE;

examples/prom/poisson_global_rom.cpp

@@ -25,7 +25,6 @@
 //
 // Online phase:  poisson_global_rom -online -f 1.15
 
-
 #include "mfem.hpp"
 #include <fstream>
 #include <iostream>

lib/algo/AdaptiveDMD.cpp

@@ -21,7 +21,8 @@ namespace CAROM {
 
 AdaptiveDMD::AdaptiveDMD(int dim, double desired_dt, std::string rbf,
                          std::string interp_method,
-                         double closest_rbf_val) : DMD(dim, desired_dt)
+                         double closest_rbf_val,
+                         Vector* state_offset) : DMD(dim, state_offset)
 {
     CAROM_VERIFY(rbf == "G" || rbf == "IQ" || rbf == "IMQ");
     CAROM_VERIFY(interp_method == "LS" || interp_method == "IDW"
@@ -33,25 +34,25 @@ AdaptiveDMD::AdaptiveDMD(int dim, double desired_dt, std::string rbf,
     d_closest_rbf_val = closest_rbf_val;
 }
 
-void AdaptiveDMD::train(double energy_fraction)
+void AdaptiveDMD::train(double energy_fraction, const Matrix* W0, double linearity_tol)
 {
     const Matrix* f_snapshots = getInterpolatedSnapshots();
     CAROM_VERIFY(f_snapshots->numColumns() > 1);
     CAROM_VERIFY(energy_fraction > 0 && energy_fraction <= 1);
     d_energy_fraction = energy_fraction;
-    constructDMD(f_snapshots, d_rank, d_num_procs);
+    constructDMD(f_snapshots, d_rank, d_num_procs, W0, linearity_tol);
 
     delete f_snapshots;
 }
 
-void AdaptiveDMD::train(int k)
+void AdaptiveDMD::train(int k, const Matrix* W0, double linearity_tol)
 {
     const Matrix* f_snapshots = getInterpolatedSnapshots();
     CAROM_VERIFY(f_snapshots->numColumns() > 1);
     CAROM_VERIFY(k > 0 && k <= f_snapshots->numColumns() - 1);
     d_energy_fraction = -1.0;
     d_k = k;
-    constructDMD(f_snapshots, d_rank, d_num_procs);
+    constructDMD(f_snapshots, d_rank, d_num_procs, W0, linearity_tol);
 
     delete f_snapshots;
 }

lib/algo/AdaptiveDMD.h

@@ -46,20 +46,28 @@ class AdaptiveDMD : public DMD
      *                            "IDW" == inverse distance weighting, "LP" == lagrangian polynomials)
      * @param[in] closest_rbf_val The RBF parameter determines the width of influence.
      *                            Set the RBF value of the nearest two parameter points to a value between 0.0 to 1.0
+     * @param[in] state_offset    The state offset.
      */
     AdaptiveDMD(int dim, double desired_dt = -1.0, std::string rbf = "G",
                 std::string interp_method = "LS",
-                double closest_rbf_val = 0.9);
+                double closest_rbf_val = 0.9,
+                Vector* state_offset = NULL);
 
     /**
      * @param[in] energy_fraction The energy fraction to keep after doing SVD.
+     * @param[in] W0              The initial basis to prepend to W.
+     * @param[in] linearity_tol   The tolerance for determining whether a column
+                                  of W is linearly independent with W0.
      */
-    void train(double energy_fraction);
+    void train(double energy_fraction, const Matrix* W0 = NULL, double linearity_tol = 0.0);
 
     /**
      * @param[in] k The number of modes (eigenvalues) to keep after doing SVD.
+     * @param[in] W0              The initial basis to prepend to W.
+     * @param[in] linearity_tol   The tolerance for determining whether a column
+                                  of W is linearly independent with W0.
      */
-    void train(int k);
+    void train(int k, const Matrix* W0 = NULL, double linearity_tol = 0.0);
 
     /**
      * @brief Get the true dt between interpolated snapshots.