Various ROM projections for steady NS toy model (#8)

* ns_rom: sv and linear operator saving.

* ns_rom: project mode to compare projection error.

* ns_rom: project - visualizing variables.

* sample generation will stop if the solution did not converge.

* sample generation: random sample generator will try another sample if not converged.

* SolveROM - enable restarting from a FOM solution projection.

* initial loading of CGOptimizer.

* MFEMROMHandler::NonlinearSolve employs cg optimizer.

* CGOptimizer Brent bug fix.

* multi_component steady ns regression test.

* sketches/ns_rom: sampling performs SVD for solution, velocity, and pressure.

* u/p snapshot generator filenames.

* change switch to if statements.

* sketches/ns_rom: TENSOR3 - vel only rom.

* sketches/ns_rom: TENSOR4 - supremizer-enriched velocity basis.

* sketches/ns_rom: TENSOR5 - petrov-galerkin, supremizer-enriched velocity projection.

include/linalg_utils.hpp

@@ -44,6 +44,9 @@ void PrintVector(const CAROM::Vector &vec,
 
 namespace mfem
 {
+
+void GramSchmidt(DenseMatrix& mat);
+
 // Compute Rt * A * P
 void RtAP(DenseMatrix& R,
         const Operator& A,
@@ -93,6 +96,42 @@ void TensorAddMultTranspose(const DenseTensor &tensor, const Vector &x, const in
 // axis 1: M_{ki} += w * T_{ijk} * x_j
 void TensorAddScaledMultTranspose(const DenseTensor &tensor, const double w, const Vector &x, const int axis, DenseMatrix &M);
 
+class CGOptimizer : public OptimizationSolver
+{
+protected:
+   const double golden_ratio = 0.5 * (1.0 + sqrt(5.0));
+   const double Cr = 1.0 - 1. / golden_ratio;
+   const double ib = 0.1;
+   const double eps = 1.0e-14;
+
+   const int N_mnbrak = 1e4;
+   const int N_para = 1e4;
+
+public:
+   CGOptimizer() : OptimizationSolver() {}
+   virtual ~CGOptimizer() {}
+
+   virtual void SetOperator(const Operator &op) override
+   {
+      oper = &op;
+      height = op.Height();
+      width = op.Width();
+
+      resvec.SetSize(height);
+   }
+
+   virtual void Mult(const Vector &rhs, Vector &sol) const;
+
+private:
+   mutable Vector resvec;
+   mutable Vector sol1, g, xi, xi1, h;
+
+   double Objective(const Vector &rhs, const Vector &sol) const;
+   void Gradient(const Vector &rhs, const Vector &sol, Vector &grad) const;
+
+   double Brent(const Vector &rhs, const Vector &xi, Vector &sol, double b0 = 1e-1, double lin_tol = 1e-2) const;
+};
+
 }
 
 #endif

include/multiblock_solver.hpp

@@ -195,7 +195,7 @@ friend class ParameterizedProblem;
 
    void AssembleROM();
 
-   virtual void Solve() = 0;
+   virtual bool Solve() = 0;
 
    virtual void InitVisualization(const std::string& output_dir = "");
    virtual void InitUnifiedParaview(const std::string &file_prefix);

include/poisson_solver.hpp

@@ -85,7 +85,7 @@ friend class ParameterizedProblem;
    virtual void BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp);
    virtual void BuildInterfaceROMElement(Array<FiniteElementSpace *> &fes_comp);
 
-   virtual void Solve();
+   virtual bool Solve();
 
    virtual void ProjectOperatorOnReducedBasis();
 

include/random_sample_generator.hpp

@@ -16,6 +16,8 @@ class RandomSampleGenerator : public SampleGenerator
 
    virtual ~RandomSampleGenerator() {}
 
+   virtual SampleGeneratorType GetType() override { return RANDOM; }
+
    // RandomSampleGenerator has the same sampling size for all parameters, equal to total samples.
    // const Array<int> GetSampleSizes() { return sampling_sizes; }
 

include/sample_generator.hpp

@@ -12,6 +12,13 @@
 
 using namespace mfem;
 
+enum SampleGeneratorType
+{
+   BASE,
+   RANDOM,
+   NUM_SAMPLE_GEN_TYPE
+};
+
 class SampleGenerator
 {
 protected:
@@ -51,6 +58,8 @@ class SampleGenerator
 
    virtual ~SampleGenerator();
 
+   virtual SampleGeneratorType GetType() { return BASE; }
+
    const int GetNumSampleParams() { return num_sampling_params; }
    const Array<int> GetSampleSizes() { return sampling_sizes; }
    const int GetTotalSampleSize() { return total_samples; }

include/steady_ns_solver.hpp

@@ -131,7 +131,7 @@ friend class SteadyNSOperator;
    virtual void SaveCompBdrROMElement(hid_t &file_id) override;
    virtual void LoadCompBdrROMElement(hid_t &file_id) override;
 
-   virtual void Solve();
+   virtual bool Solve();
 
    virtual void ProjectOperatorOnReducedBasis();
 

include/stokes_solver.hpp

@@ -185,7 +185,7 @@ friend class ParameterizedProblem;
    virtual void BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp);
    virtual void BuildInterfaceROMElement(Array<FiniteElementSpace *> &fes_comp);
 
-   virtual void Solve();
+   virtual bool Solve();
    virtual void Solve_obsolete();
 
    virtual void ProjectOperatorOnReducedBasis();