AdvDiffSolver flow-solver upgrade + alternating-Schwarz ROM for SteadyNS

CMakeLists.txt

@@ -211,6 +211,7 @@ set(scaleupROMObj_SOURCES
 
   include/steady_ns_solver.hpp
   src/steady_ns_solver.cpp
+  src/steady_ns_solver_schwarz.cpp
 
   include/advdiff_solver.hpp
   src/advdiff_solver.cpp

bin/main.cpp

@@ -34,6 +34,7 @@ int main(int argc, char *argv[])
       else if (mode == "aux_train_rom")   AuxiliaryTrainROM(MPI_COMM_WORLD);
       else if (mode == "train_eqp")       TrainEQP(MPI_COMM_WORLD);
       else if (mode == "single_run")      SingleRun(MPI_COMM_WORLD, output_file);
+      else if (mode == "schwarz_rom")     SingleSchwarzRun(MPI_COMM_WORLD, output_file);
       else if (mode == "print_eqp")       PrintEQPCoords(MPI_COMM_WORLD);
       else
       {

include/advdiff_solver.hpp

@@ -6,7 +6,7 @@
 #define SCALEUPROM_ADVDIFF_SOLVER_HPP
 
 #include "poisson_solver.hpp"
-#include "stokes_solver.hpp"
+#include "steady_ns_solver.hpp"
 
 // By convention we only use mfem namespace as default, not CAROM.
 using namespace mfem;
@@ -26,7 +26,9 @@ friend class ParameterizedProblem;
    /*
       flow solver to obtain the prescribed velocity field. both StokesSolver / SteadyNSSolver can be used.
    */
-   StokesSolver *stokes_solver = NULL;
+   std::string flow_solver_type = "";
+   int flow_solver_order = -1;
+   StokesSolver *flow_solver = NULL;
    bool load_flow = false;
    bool save_flow = false;
    std::string flow_file = "";
@@ -39,7 +41,7 @@ friend class ParameterizedProblem;
    Array<GridFunction *> global_flow_visual;
 
 public:
-   AdvDiffSolver();
+   AdvDiffSolver(TopologyHandler *input_topol_handler=NULL);
 
    virtual ~AdvDiffSolver();
 
@@ -52,15 +54,15 @@ friend class ParameterizedProblem;
 
    void SetFlowAtSubdomain(std::function<void(const Vector &, double, Vector &)> F, const int m=-1);
 
-   void SetParameterizedProblem(ParameterizedProblem *problem) override;
+   bool SetParameterizedProblem(ParameterizedProblem *problem) override;
 
    void SaveVisualization() override;
 
 protected:
    void SetMUMPSSolver() override;
 
 private:
-   void GetFlowField(ParameterizedProblem *flow_problem);
+   bool GetFlowField(ParameterizedProblem *flow_problem);
 };
 
 #endif

include/component_topology_handler.hpp

@@ -83,6 +83,13 @@ class ComponentTopologyHandler : public TopologyHandler
 
       std::unordered_map<int,int> vtx2to1;             // vertex mapping from component 2 to component 1.
       Array2D<int> be_pairs;  // boundary element pairs between component 1 and 2.
+
+      PortData() {}
+
+      PortData(const PortData &input)
+         : Component1(input.Component1), Component2(input.Component2),
+           Attr1(input.Attr1), Attr2(input.Attr2), vtx2to1(input.vtx2to1),
+           be_pairs(input.be_pairs) {}
    };
 
 protected:
@@ -133,6 +140,19 @@ class ComponentTopologyHandler : public TopologyHandler
 public:
    ComponentTopologyHandler();
 
+   // Get a subset of the mesh array as a new ComponentTopologyHandler.
+   // Assumes meshes are arranged in an N x N square array.
+   // Extracts an M x M subset starting at (i0, j0).
+   // Internal ports within the subset are included.
+   // External ports connecting to meshes outside the subset are excluded;
+   // their port attributes (formerly interior to the global domain) are
+   // registered into bdr_attributes as new global boundary attributes of the subset mesh.
+   // Subset mesh boundaries that lay on the global boundary inherit the
+   // corresponding global boundary attributes.
+   // Stores subset mesh to original mesh mapping.
+   ComponentTopologyHandler(ComponentTopologyHandler* global, const int i0, const int j0,
+                            const int N, const int M, Array<int> &subset_to_orig);
+
    virtual ~ComponentTopologyHandler();
 
    // access

include/linelast_solver.hpp

@@ -57,7 +57,7 @@ class LinElastSolver : public MultiBlockSolver
    VectorCoefficient *init_x = NULL;
 
 public:
-   LinElastSolver();
+   LinElastSolver(TopologyHandler *input_topol_handler=NULL);
 
    virtual ~LinElastSolver();
 
@@ -86,6 +86,7 @@ class LinElastSolver : public MultiBlockSolver
    virtual void SetupBCVariables() override;
    virtual void SetupIC(std::function<void(const Vector &, double, Vector &)> F);
    virtual void AddBCFunction(std::function<void(const Vector &, double, Vector &)> F, const int battr = -1);
+   virtual void AddBCFunction(const Vector &F, const int battr = -1);
    virtual void AddRHSFunction(std::function<void(const Vector &, double, Vector &)> F);
    virtual bool BCExistsOnBdr(const int &global_battr_idx);
    virtual void SetupBCOperators();
@@ -101,7 +102,7 @@ class LinElastSolver : public MultiBlockSolver
 
    void SanityCheckOnCoeffs();
 
-   virtual void SetParameterizedProblem(ParameterizedProblem *problem);
+   virtual bool SetParameterizedProblem(ParameterizedProblem *problem);
 };
 
 #endif

include/main_workflow.hpp

@@ -32,6 +32,9 @@ void FindSnapshotFilesForBasis(const BasisTag &basis_tag, const std::string &def
 // return relative error if comparing solution.
 double SingleRun(MPI_Comm comm, const std::string output_file = "");
 
+// return relative error if comparing solution.
+double SingleSchwarzRun(MPI_Comm comm, const std::string &output_file = "");
+
 // Auxiliary function to print out EQP point coordinates
 void PrintEQPCoords(MPI_Comm comm);
 

include/multiblock_solver.hpp

@@ -116,7 +116,7 @@ friend class ParameterizedProblem;
    ROMLinearElement *rom_elems = NULL;
 
 public:
-   MultiBlockSolver();
+   MultiBlockSolver(TopologyHandler *input_topol_handler=NULL);
 
    virtual ~MultiBlockSolver();
 
@@ -254,14 +254,25 @@ friend class ParameterizedProblem;
    virtual void SaveBasisVisualization()
    { rom_handler->SaveBasisVisualization(fes, var_names); }
 
-   virtual void SetParameterizedProblem(ParameterizedProblem *problem);
+   virtual bool SetParameterizedProblem(ParameterizedProblem *problem);
 
    void ComputeSubdomainErrorAndNorm(GridFunction *fom_sol, GridFunction *rom_sol, double &error, double &norm);
    void ComputeRelativeError(Array<GridFunction *> fom_sols, Array<GridFunction *> rom_sols, Vector &error);
    void CompareSolution(BlockVector &test_U, Vector &error);
 
    virtual void SaveEQPCoords(const std::string &filename) {}
 
+   bool IsBdrTypeDefined()
+   {
+      for (int k = 0; k < bdr_type.Size(); k++)
+         if (bdr_type[k] == BoundaryType::NUM_BDR_TYPE)
+            return false;
+
+      return true;
+   }
+
+   void PrintConfiguration() const;
+
 protected:
    virtual void AssembleROMMat(BlockMatrix &romMat);
 };

include/parameterized_problem.hpp

@@ -218,6 +218,8 @@ enum class BoundaryType
    NUM_BDR_TYPE
 };
 
+const std::string PrintBoundaryType(const BoundaryType bdr_type);
+
 class ParameterizedProblem
 {
 friend class SampleGenerator;

include/poisson_solver.hpp

@@ -50,7 +50,7 @@ friend class ParameterizedProblem;
    double kappa = -1.0;
 
 public:
-   PoissonSolver();
+   PoissonSolver(TopologyHandler *input_topol_handler=NULL);
 
    virtual ~PoissonSolver();
 
@@ -98,7 +98,7 @@ friend class ParameterizedProblem;
 
    void SanityCheckOnCoeffs();
 
-   virtual void SetParameterizedProblem(ParameterizedProblem *problem);
+   virtual bool SetParameterizedProblem(ParameterizedProblem *problem);
 
 protected:
    virtual void SetMUMPSSolver();

include/steady_ns_solver.hpp

@@ -99,6 +99,8 @@ class SteadyNSTensorROM : public SteadyNSROM
 
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual Operator &GetGradient(const Vector &x) const;
+
+   void SaveOperator(std::string filename);
 };
 
 class SteadyNSEQPROM : public SteadyNSROM
@@ -166,7 +168,7 @@ friend class SteadyNSOperator;
    NewtonSolver *newton_solver = NULL;
 
 public:
-   SteadyNSSolver();
+   SteadyNSSolver(TopologyHandler *input_topol_handler=NULL);
 
    virtual ~SteadyNSSolver();
 
@@ -198,6 +200,12 @@ friend class SteadyNSOperator;
 
    void SaveEQPCoords(const std::string &filename) override;
 
+   void SchwarzROM(const int M, const int N, ParameterizedProblem *problem,
+                   double &solve_time, int &num_solve, Array<double> &error_hist,
+                   const int maxIter, const double threshold=1e-15,
+                   const int plateau_track=3, const double plateau_range=1e-1,
+                   const bool use_restart=false, const double initial_tol=-1.0);
+
 private:
    DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace);
 
@@ -210,6 +218,10 @@ friend class SteadyNSOperator;
    void SaveEQPElems(const std::string &filename);
    void LoadEQPElems(const std::string &filename);
    void AssembleROMEQPOper();
+
+protected:
+   void SetupSubsetRHSBCOperators(const Array<int> &meshes, const Array<int> &battrs,
+                                  const Array<VectorGridFunctionCoefficient *> &bfuncs);
 };
 
 #endif

include/stokes_solver.hpp

@@ -134,7 +134,7 @@ friend class ParameterizedProblem;
    double kappa = -1.0;
 
 public:
-   StokesSolver();
+   StokesSolver(TopologyHandler *input_topol_handler=NULL);
 
    virtual ~StokesSolver();
 
@@ -197,7 +197,7 @@ friend class ParameterizedProblem;
 
    void SanityCheckOnCoeffs();
 
-   virtual void SetParameterizedProblem(ParameterizedProblem *problem) override;
+   virtual bool SetParameterizedProblem(ParameterizedProblem *problem) override;
 
    // to ensure incompressibility for the problems with all velocity dirichlet bc.
    void SetComplementaryFlux(const Array<bool> nz_dbcs);
@@ -216,6 +216,19 @@ friend class ParameterizedProblem;
    double ComputeBEIntegral(const FiniteElement &el, ElementTransformation &Tr, Coefficient &Q);
    void ComputeBEIntegral(const FiniteElement &el, ElementTransformation &Tr,
                            VectorCoefficient &Q, Vector &result);
+
+protected:
+   // For an M x M subset (starting at (i0, j0)) extracted from the assumed
+   // N x N square mesh array: when all subset boundary conditions are velocity
+   // Dirichlet, the prescribed velocity must satisfy the compatibility
+   // condition int_{dOmega} u . n dA = 0 for the incompressible Stokes problem
+   // to be well-posed. This routine computes the net flux from the prescribed
+   // BC and adds a complementary flux on the subset boundary to enforce that
+   // constraint.
+   void SetSubsetComplementaryFlux(const int N, const int M, const int i0, const int j0,
+                                   const Array<int> &subset_bdr_attributes,
+                                   const Array<BoundaryType> &subset_bdrtype,
+                                   const ParameterizedProblem *problem);
 };
 
 #endif

include/unsteady_ns_solver.hpp

@@ -74,7 +74,7 @@ friend class SteadyNSOperator;
    double times[10];
 
 public:
-   UnsteadyNSSolver();
+   UnsteadyNSSolver(TopologyHandler *input_topol_handler=NULL);
 
    virtual ~UnsteadyNSSolver();
 
@@ -93,7 +93,7 @@ friend class SteadyNSOperator;
    using MultiBlockSolver::SaveVisualization;
    void SaveVisualization(const int step, const double time) override;
 
-   void SetParameterizedProblem(ParameterizedProblem *problem) override;
+   bool SetParameterizedProblem(ParameterizedProblem *problem) override;
 
    BlockVector* PrepareSnapshots(std::vector<BasisTag> &basis_tags) override;
 

sketches/CMakeLists.txt

@@ -1,4 +1,5 @@
 add_executable(poisson poisson.cpp $<TARGET_OBJECTS:scaleupROMObj>)
+add_executable(phyml_poisson phyml_poisson.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(poisson_serial poisson_serial.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(poisson_1d poisson_1d.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(multidomain multidomain.cpp $<TARGET_OBJECTS:scaleupROMObj>)