Unsteady Navier-Stokes solver (#45)

* steady ns solver with LF integrator, interface not yet tested.

* steady ns solver with LF integrator verified.

* initial loading for UnsteadyNSSolver

* initial loading of unsteady ns time integration. not tested yet

* UnsteadyNSSolver mms verified.

* time integration input

* MultiBlockSolver visualization option consolidated.

* time interval for visualization

* name changes for flow parameterized problems.

* backward facing step flow problem

* visualization, compute cfl, sanity check.

* SteadyNSSolver: LF integrator should be added for neumann condition as well.

* UnsteadyNSSolver::SetTime

* report interval, fix in cfl.

* general functions now have time as input.

* added restart I/O feature

* periodic flow past array problem

.github/workflows/ci.yml

@@ -77,6 +77,10 @@ jobs:
         run: |
             cd ${GITHUB_WORKSPACE}/build/test
             ./steady_ns_dd_mms
+      - name: Test UnsteadyNS DD solver
+        run: |
+            cd ${GITHUB_WORKSPACE}/build/test
+            ./unsteady_ns_dd_mms
       - name: Test LinearElastic DD solver
         run: |
             cd ${GITHUB_WORKSPACE}/build/test

CMakeLists.txt

@@ -160,6 +160,8 @@ set(scaleupROMObj_SOURCES
 
   # include/navier_solver.hpp
   # src/navier_solver.cpp
+  include/unsteady_ns_solver.hpp
+  src/unsteady_ns_solver.cpp
 
   include/topology_handler.hpp
   src/topology_handler.cpp

include/advdiff_solver.hpp

@@ -50,7 +50,7 @@ friend class ParameterizedProblem;
 
    bool Solve() override;
 
-   void SetFlowAtSubdomain(std::function<void(const Vector &, Vector &)> F, const int m=-1);
+   void SetFlowAtSubdomain(std::function<void(const Vector &, double, Vector &)> F, const int m=-1);
 
    void SetParameterizedProblem(ParameterizedProblem *problem) override;
 

include/etc.hpp

@@ -27,4 +27,15 @@ inline void DeletePointers(Array2D<T*> &ptr_array)
 
 bool FileExists(const std::string& name);
 
+template<typename ... Args>
+std::string string_format( const std::string& format, Args ... args )
+{
+    int size_s = std::snprintf( nullptr, 0, format.c_str(), args ... ) + 1; // Extra space for '\0'
+    if( size_s <= 0 ){ throw std::runtime_error( "Error during formatting." ); }
+    auto size = static_cast<size_t>( size_s );
+    std::unique_ptr<char[]> buf( new char[ size ] );
+    std::snprintf( buf.get(), size, format.c_str(), args ... );
+    return std::string( buf.get(), buf.get() + size - 1 ); // We don't want the '\0' inside
+}
+
 #endif

include/linelast_solver.hpp

@@ -40,10 +40,10 @@ class LinElastSolver : public MultiBlockSolver
    Array<BilinearForm *> as;
 
    // Lame constants for each subdomain, global boundary attribute ordering
-   Array<ConstantCoefficient *> lambda_c;
-   Array<ConstantCoefficient *> mu_c;
-   Array<VectorFunctionCoefficient *> bdr_coeffs;
-   Array<VectorFunctionCoefficient *> rhs_coeffs;
+   Array<Coefficient *> lambda_c;
+   Array<Coefficient *> mu_c;
+   Array<VectorCoefficient *> bdr_coeffs;
+   Array<VectorCoefficient *> rhs_coeffs;
 
    // DG parameters specific to linear elasticity equation.
    double alpha = -1.0;
@@ -54,7 +54,7 @@ class LinElastSolver : public MultiBlockSolver
    double mu = 1.0;
 
    // Initial positions
-   VectorFunctionCoefficient *init_x = NULL;
+   VectorCoefficient *init_x = NULL;
 
 public:
    LinElastSolver();
@@ -84,9 +84,9 @@ class LinElastSolver : public MultiBlockSolver
    virtual bool Solve();
 
    virtual void SetupBCVariables() override;
-   virtual void SetupIC(std::function<void(const Vector &, Vector &)> F);
-   virtual void AddBCFunction(std::function<void(const Vector &, Vector &)> F, const int battr = -1);
-   virtual void AddRHSFunction(std::function<void(const Vector &, Vector &)> F);
+   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 AddRHSFunction(std::function<void(const Vector &, double, Vector &)> F);
    virtual bool BCExistsOnBdr(const int &global_battr_idx);
    virtual void SetupBCOperators();
    virtual void SetupRHSBCOperators();

include/main_workflow.hpp

@@ -9,8 +9,8 @@
 #include "multiblock_solver.hpp"
 #include "random_sample_generator.hpp"
 
-double dbc2(const Vector &);
-double dbc4(const Vector &);
+double dbc2(const Vector &, double t);
+double dbc4(const Vector &, double t);
 void RunExample();
 
 MultiBlockSolver* InitSolver();

include/multiblock_solver.hpp

@@ -79,18 +79,29 @@ friend class ParameterizedProblem;
    std::string sol_dir = ".";
    std::string sol_prefix;
 
+   /* visualization options */
+   struct VisualizationOption {
+      bool save = false;
+      bool unified_view = false;
+      std::string dir = ".";
+      std::string prefix;
+      /* for individual_paraview, support partial visualization. */
+      int domain_offset = -1;
+      int domain_interval = -1;
+
+      /* time-dependent options */
+      int time_interval = -1;
+
+      /* visualizing rom solution error */
+      bool save_error = false;
+   } visual;
+
    // visualization variables
-   bool save_visual = false;
-   bool unified_paraview = false;
-   int visual_offset = -1, visual_freq = -1;    // for individual_paraview, support partial visualization.
-   std::string visual_dir = ".";
-   std::string visual_prefix;
    Array<ParaViewDataCollection *> paraviewColls;
    // Used only for the unified visualization. Size(num_var).
    Array<FiniteElementSpace *> global_fes;
    Array<GridFunction *> global_us_visual;
 
-   bool visualize_error = false;
    Array<GridFunction *> error_visual;    // Size(num_var * numSub)
    Array<GridFunction *> global_error_visual;   // point-wise error visualization
 
@@ -130,9 +141,9 @@ friend class ParameterizedProblem;
    ROMHandlerBase* GetROMHandler() const { return rom_handler; }
    TopologyHandler* GetTopologyHandler() const { return topol_handler; }
    const TrainMode GetTrainMode() { return train_mode; }
-   const bool IsVisualizationSaved() const { return save_visual; }
+   const bool IsVisualizationSaved() const { return visual.save; }
    const std::string GetSolutionFilePrefix() const { return sol_prefix; }
-   const std::string GetVisualizationPrefix() const { return visual_prefix; }
+   const std::string GetVisualizationPrefix() const { return visual.prefix; }
    const TopologyHandlerMode GetTopologyMode() const { return topol_mode; }
    ParaViewDataCollection* GetParaViewColl(const int &k) { return paraviewColls[k]; }
    BlockVector* GetSolution() { return U; }
@@ -147,11 +158,11 @@ friend class ParameterizedProblem;
    void SortByVariables(BlockVector &by_sub, BlockVector &by_var);
 
    virtual void SetupBCVariables();
-   virtual void AddBCFunction(std::function<double(const Vector &)> F, const int battr = -1)
+   virtual void AddBCFunction(std::function<double(const Vector &, double)> F, const int battr = -1)
    { mfem_error("Abstract method MultiBlockSolver::AddBCFunction!\n"); }
    virtual void AddBCFunction(const double &F, const int battr = -1)
    { mfem_error("Abstract method MultiBlockSolver::AddBCFunction!\n"); }
-   virtual void AddBCFunction(std::function<void(const Vector &, Vector &)> F, const int battr = -1)
+   virtual void AddBCFunction(std::function<void(const Vector &, double, Vector &)> F, const int battr = -1)
    { mfem_error("Abstract method MultiBlockSolver::AddBCFunction!\n"); }
    virtual void AddBCFunction(const Vector &F, const int battr = -1)
    { mfem_error("Abstract method MultiBlockSolver::AddBCFunction!\n"); }
@@ -220,10 +231,15 @@ friend class ParameterizedProblem;
    virtual void InitVisualization(const std::string& output_dir = "");
    virtual void InitUnifiedParaview(const std::string &file_prefix);
    virtual void InitIndividualParaview(const std::string &file_prefix);
+   /* time-independent visualization */
    virtual void SaveVisualization();
+   /* time-dependent visualization */
+   virtual void SaveVisualization(const int step, const double time);
 
    void SaveSolution(std::string filename = "");
+   void SaveSolutionWithTime(const std::string filename, const int step, const double time);
    void LoadSolution(const std::string &filename);
+   void LoadSolutionWithTime(const std::string &filename, int &step, double &time);
    void CopySolution(BlockVector *input_sol);
 
    virtual void AllocateROMTensorElems()