Advdiff update

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,8 @@ 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 = "";
+   StokesSolver *flow_solver = NULL;
    bool load_flow = false;
    bool save_flow = false;
    std::string flow_file = "";
@@ -52,15 +53,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/linelast_solver.hpp

@@ -101,7 +101,7 @@ class LinElastSolver : public MultiBlockSolver
 
    void SanityCheckOnCoeffs();
 
-   virtual void SetParameterizedProblem(ParameterizedProblem *problem);
+   virtual bool SetParameterizedProblem(ParameterizedProblem *problem);
 };
 
 #endif

include/multiblock_solver.hpp

@@ -254,7 +254,7 @@ 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);

include/poisson_solver.hpp

@@ -98,7 +98,7 @@ friend class ParameterizedProblem;
 
    void SanityCheckOnCoeffs();
 
-   virtual void SetParameterizedProblem(ParameterizedProblem *problem);
+   virtual bool SetParameterizedProblem(ParameterizedProblem *problem);
 
 protected:
    virtual void SetMUMPSSolver();

include/stokes_solver.hpp

@@ -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);

include/unsteady_ns_solver.hpp

@@ -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;
 

src/advdiff_solver.cpp

@@ -25,15 +25,17 @@ AdvDiffSolver::AdvDiffSolver()
    load_flow = config.GetOption<bool>("adv-diff/load_flow", false);
    if (save_flow || load_flow)
       flow_file = config.GetRequiredOption<std::string>("adv-diff/flow_file");
+
+   flow_solver_type = config.GetOption<std::string>("adv-diff/flow_solver_type", "stokes");
 }
 
 AdvDiffSolver::~AdvDiffSolver()
 {
    DeletePointers(flow_coeffs);
-   if (!stokes_solver) DeletePointers(flow_visual);
+   if (!flow_solver) DeletePointers(flow_visual);
    DeletePointers(global_flow_visual);
    delete global_flow_fes;
-   delete stokes_solver;
+   delete flow_solver;
 }
 
 void AdvDiffSolver::BuildDomainOperators()
@@ -169,12 +171,14 @@ void AdvDiffSolver::SetFlowAtSubdomain(std::function<void(const Vector &, double
       flow_coeffs[m] = new VectorFunctionCoefficient(dim, F);
 }
 
-void AdvDiffSolver::SetParameterizedProblem(ParameterizedProblem *problem)
+bool AdvDiffSolver::SetParameterizedProblem(ParameterizedProblem *problem)
 {
+   bool success = true;
    if (!function_factory::advdiff_problem::analytic_flow)
-      GetFlowField(function_factory::advdiff_problem::flow_problem);
+      success = GetFlowField(function_factory::advdiff_problem::flow_problem);
 
-   PoissonSolver::SetParameterizedProblem(problem);
+   success = success && (PoissonSolver::SetParameterizedProblem(problem));
+   return success;
 }
 
 void AdvDiffSolver::SaveVisualization()
@@ -188,8 +192,8 @@ void AdvDiffSolver::SaveVisualization()
    flow_visual.SetSize(numSub);
    flow_visual = NULL;
 
-   const int stokes_numvar = (stokes_solver) ? stokes_solver->GetNumVar() : 0;
-   if (!stokes_solver)
+   const int stokes_numvar = (flow_solver) ? flow_solver->GetNumVar() : 0;
+   if (!flow_solver)
    {
       flow_fes.SetSize(numSub);
       for (int m = 0; m < numSub; m++)
@@ -201,8 +205,8 @@ void AdvDiffSolver::SaveVisualization()
 
    for (int m = 0; m < numSub; m++)
    {
-      if (stokes_solver)
-         flow_visual[m] = stokes_solver->GetGridFunction(m * stokes_numvar);
+      if (flow_solver)
+         flow_visual[m] = flow_solver->GetGridFunction(m * stokes_numvar);
       else
       {
          assert(flow_coeffs[m]);
@@ -235,43 +239,51 @@ void AdvDiffSolver::SetMUMPSSolver()
    mumps->SetOperator(*globalMat_hypre);
 }
 
-void AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem)
+bool AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem)
 {
    assert(flow_problem);
    mfem_warning("AdvDiffSolver: Obtaining flow field. This may take a while depending on the domain size.\n");
 
-   stokes_solver = new StokesSolver;
-   stokes_solver->InitVariables();
-   if (use_rom) stokes_solver->InitROMHandler();
-   stokes_solver->SetSolutionSaveMode(save_flow);
+   if (flow_solver_type == "stokes")
+      flow_solver = new StokesSolver;
+   else if (flow_solver_type == "steady-ns")
+      flow_solver = new SteadyNSSolver;
+   else
+      mfem_error("AdvDiffSolver::GetFlowField - Unknown flow solver type!\n");
+
+   flow_solver->InitVariables();
+   if (use_rom) flow_solver->InitROMHandler();
+   flow_solver->SetSolutionSaveMode(save_flow);
 
    bool flow_loaded = false;
    if (load_flow && FileExists(flow_file))
    {
-      stokes_solver->LoadSolution(flow_file);
+      flow_solver->LoadSolution(flow_file);
       flow_loaded = true;
    }
    else
    {
-      stokes_solver->SetParameterizedProblem(flow_problem);
+      flow_solver->SetParameterizedProblem(flow_problem);
       // currently only support FOM.
-      stokes_solver->BuildOperators();
-      stokes_solver->SetupBCOperators();
-      stokes_solver->Assemble();
-      stokes_solver->Solve();
+      flow_solver->BuildOperators();
+      flow_solver->SetupBCOperators();
+      flow_solver->Assemble();
+      flow_loaded = flow_solver->Solve();
    }
 
    if (save_flow && (!flow_loaded))
-      stokes_solver->SaveSolution(flow_file);
+      flow_solver->SaveSolution(flow_file);
 
    DeletePointers(flow_coeffs);
-   const int stokes_numvar = stokes_solver->GetNumVar();
+   const int stokes_numvar = flow_solver->GetNumVar();
    for (int m = 0; m < numSub; m++)
-      flow_coeffs[m] = new VectorGridFunctionCoefficient(stokes_solver->GetGridFunction(m * stokes_numvar));
+      flow_coeffs[m] = new VectorGridFunctionCoefficient(flow_solver->GetGridFunction(m * stokes_numvar));
 
    /*
       VectorGridFunctionCoefficient does not own the grid function,
       and it requires the grid function for its lifetime.
-      Thus stokes_solver will be deleted at ~AdvDiffSolver().
+      Thus flow_solver will be deleted at ~AdvDiffSolver().
    */
+
+   return flow_loaded;
 }

src/linelast_solver.cpp

@@ -440,7 +440,7 @@ void LinElastSolver::SetupDomainBCOperators()
    }
 }
 
-void LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem)
+bool LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem)
 {
    /* set up boundary types */
    MultiBlockSolver::SetParameterizedProblem(problem);
@@ -483,6 +483,8 @@ void LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem)
    {
       SetupIC(*(problem->general_vector_ptr[0]));
    }
+   // LinElastSolver does not fail in SetParameterizedProblem.
+   return true;
 }
 
 void LinElastSolver::ProjectOperatorOnReducedBasis()

src/main_workflow.cpp

@@ -165,7 +165,22 @@ void GenerateSamples(MPI_Comm comm)
          test->InitROMHandler();
 
       problem->SetSingleRun();
-      test->SetParameterizedProblem(problem);
+      bool param_set = test->SetParameterizedProblem(problem);
+      if (!param_set)
+      {
+         if (sample_generator->GetType() != SampleGeneratorType::RANDOM)
+            mfem_error("GenerateSamples failed at SetParameterizedProblem, and parameter values are fixed!\n");
+         delete test;
+
+         sample_generator->SetSampleParams(s);
+         test = InitSolver();
+         test->InitVariables();
+         if (test->UseRom())
+            test->InitROMHandler();
+
+         problem->SetSingleRun();
+         param_set = test->SetParameterizedProblem(problem);
+      }
 
       int file_idx = s + sample_generator->GetFileOffset();
       const std::string visual_path = sample_generator->GetSamplePath(file_idx, test->GetVisualizationPrefix());
@@ -427,7 +442,7 @@ void BuildROM(MPI_Comm comm)
 
    // The ROM operator will be built based on the parameter specified for single-run.
    problem->SetSingleRun();
-   test->SetParameterizedProblem(problem);
+   assert(test->SetParameterizedProblem(problem));
 
    // TODO: there are skippable operations depending on rom/fom mode.
    test->BuildOperators();
@@ -504,7 +519,7 @@ double SingleRun(MPI_Comm comm, const std::string output_file)
    std::string solveType = (test->UseRom()) ? "ROM" : "FOM";
 
    problem->SetSingleRun();
-   test->SetParameterizedProblem(problem);
+   assert(test->SetParameterizedProblem(problem));
 
    // TODO: there are skippable operations depending on rom/fom mode.
    test->BuildRHSOperators();
@@ -721,7 +736,7 @@ void PrintEQPCoords(MPI_Comm comm)
    StopWatch solveTimer;
 
    problem->SetSingleRun();
-   test->SetParameterizedProblem(problem);
+   assert(test->SetParameterizedProblem(problem));
 
    // TODO: there are skippable operations depending on rom/fom mode.
    test->BuildRHSOperators();

src/multiblock_solver.cpp

@@ -492,7 +492,7 @@ void MultiBlockSolver::SaveVisualization(const int step, const double time)
    SaveVisualization();
 }
 
-void MultiBlockSolver::SetParameterizedProblem(ParameterizedProblem *problem)
+bool MultiBlockSolver::SetParameterizedProblem(ParameterizedProblem *problem)
 {
    assert(bdr_type.Size() == global_bdr_attributes.Size());
    for (int b = 0; b < global_bdr_attributes.Size(); b++)
@@ -502,6 +502,8 @@ void MultiBlockSolver::SetParameterizedProblem(ParameterizedProblem *problem)
 
       bdr_type[b] = problem->bdr_type[idx];
    }
+   // MultiBlockSolver does not fail in SetParameterizedProblem.
+   return true;
 }
 
 void MultiBlockSolver::SaveSolution(std::string filename)

src/poisson_solver.cpp

@@ -547,7 +547,7 @@ void PoissonSolver::SanityCheckOnCoeffs()
       MFEM_WARNING("All bc coefficients are NULL, meaning there is no Dirichlet BC. Make sure to set bc coefficients before SetupBCOperator.\n");
 }
 
-void PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem)
+bool PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem)
 {
    /* set up boundary types */
    MultiBlockSolver::SetParameterizedProblem(problem);
@@ -582,6 +582,8 @@ void PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem)
       AddRHSFunction(*(problem->scalar_rhs_ptr));
    else
       AddRHSFunction(0.0);
+   // PoissonSolver does not fail in SetParameterizedProblem.
+   return true;
 }
 
 void PoissonSolver::SetMUMPSSolver()

src/rom_handler.cpp

@@ -536,10 +536,10 @@ void MFEMROMHandler::Solve(BlockVector &rhs, BlockVector &sol)
 {
    assert(operator_loaded);
 
-   int maxIter = config.GetOption<int>("solver/max_iter", 10000);
-   double rtol = config.GetOption<double>("solver/relative_tolerance", 1.e-15);
-   double atol = config.GetOption<double>("solver/absolute_tolerance", 1.e-15);
-   int print_level = config.GetOption<int>("solver/print_level", 0);   
+   int maxIter = config.GetOption<int>("rom_solver/max_iter", 10000);
+   double rtol = config.GetOption<double>("rom_solver/relative_tolerance", 1.e-15);
+   double atol = config.GetOption<double>("rom_solver/absolute_tolerance", 1.e-15);
+   int print_level = config.GetOption<int>("rom_solver/print_level", 0);   
    std::string prec_str = config.GetOption<std::string>("model_reduction/preconditioner", "none");
 
    if (linsol_type == SolverType::DIRECT)
@@ -637,8 +637,8 @@ void MFEMROMHandler::NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec
 
    printf("Solve ROM.\n");
    reduced_sol = new BlockVector(rom_block_offsets);
-   bool use_restart = config.GetOption<bool>("solver/use_restart", false);
-   double amp = config.GetOption<double>("solver/initial_random_perturbation", 1.0e-5);
+   bool use_restart = config.GetOption<bool>("rom_solver/use_restart", false);
+   double amp = config.GetOption<double>("rom_solver/initial_random_perturbation", 1.0e-5);
    if (use_restart)
       ProjectGlobalToDomainBasis(U, reduced_sol);
    else
@@ -647,15 +647,15 @@ void MFEMROMHandler::NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec
          (*reduced_sol)(k) = amp * UniformRandom();
    }
 
-   int maxIter = config.GetOption<int>("solver/max_iter", 100);
-   double rtol = config.GetOption<double>("solver/relative_tolerance", 1.e-10);
-   double atol = config.GetOption<double>("solver/absolute_tolerance", 1.e-10);
-   int print_level = config.GetOption<int>("solver/print_level", 0);   
+   int maxIter = config.GetOption<int>("rom_solver/max_iter", 100);
+   double rtol = config.GetOption<double>("rom_solver/relative_tolerance", 1.e-10);
+   double atol = config.GetOption<double>("rom_solver/absolute_tolerance", 1.e-10);
+   int print_level = config.GetOption<int>("rom_solver/print_level", 0);   
 
-   int jac_maxIter = config.GetOption<int>("solver/jacobian/max_iter", 10000);
-   double jac_rtol = config.GetOption<double>("solver/jacobian/relative_tolerance", 1.e-10);
-   double jac_atol = config.GetOption<double>("solver/jacobian/absolute_tolerance", 1.e-10);
-   int jac_print_level = config.GetOption<int>("solver/jacobian/print_level", -1);
+   int jac_maxIter = config.GetOption<int>("rom_solver/jacobian/max_iter", 10000);
+   double jac_rtol = config.GetOption<double>("rom_solver/jacobian/relative_tolerance", 1.e-10);
+   double jac_atol = config.GetOption<double>("rom_solver/jacobian/absolute_tolerance", 1.e-10);
+   int jac_print_level = config.GetOption<int>("rom_solver/jacobian/print_level", -1);
    std::string prec_str = config.GetOption<std::string>("model_reduction/preconditioner", "none");
    if (prec_str != "none") assert(prec);
 

src/stokes_solver.cpp

@@ -1079,7 +1079,7 @@ void StokesSolver::SanityCheckOnCoeffs()
       MFEM_WARNING("All velocity bc coefficients are NULL, meaning there is no Dirichlet BC. Make sure to set bc coefficients before SetupBCOperator.\n");
 }
 
-void StokesSolver::SetParameterizedProblem(ParameterizedProblem *problem)
+bool StokesSolver::SetParameterizedProblem(ParameterizedProblem *problem)
 {
    /* set up boundary types */
    MultiBlockSolver::SetParameterizedProblem(problem);
@@ -1143,6 +1143,9 @@ void StokesSolver::SetParameterizedProblem(ParameterizedProblem *problem)
       }
       SetComplementaryFlux(nz_dbcs);
    }
+
+   // StokesSolver does not fail at SetParameterizedProblem.
+   return true;
 }
 
 BlockMatrix* StokesSolver::FormBlockMatrix(
@@ -1445,4 +1448,4 @@ void StokesSolver::ComputeBEIntegral(
       Qvec *= Tr.Weight() * ip.weight;
       result += Qvec;
    }
-}
+}