diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index 62a39533..c39691d7 100644 --- a/include/advdiff_solver.hpp +++ b/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,7 +53,7 @@ friend class ParameterizedProblem; void SetFlowAtSubdomain(std::function F, const int m=-1); - void SetParameterizedProblem(ParameterizedProblem *problem) override; + bool SetParameterizedProblem(ParameterizedProblem *problem) override; void SaveVisualization() override; @@ -60,7 +61,7 @@ friend class ParameterizedProblem; void SetMUMPSSolver() override; private: - void GetFlowField(ParameterizedProblem *flow_problem); + bool GetFlowField(ParameterizedProblem *flow_problem); }; #endif diff --git a/include/linelast_solver.hpp b/include/linelast_solver.hpp index 2ea41596..e499aae4 100644 --- a/include/linelast_solver.hpp +++ b/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 diff --git a/include/multiblock_solver.hpp b/include/multiblock_solver.hpp index ea327e0d..d7f3ca4c 100644 --- a/include/multiblock_solver.hpp +++ b/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 fom_sols, Array rom_sols, Vector &error); diff --git a/include/poisson_solver.hpp b/include/poisson_solver.hpp index 9aa72f86..115dd1ee 100644 --- a/include/poisson_solver.hpp +++ b/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(); diff --git a/include/stokes_solver.hpp b/include/stokes_solver.hpp index e1171956..614b41b9 100644 --- a/include/stokes_solver.hpp +++ b/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 nz_dbcs); diff --git a/include/unsteady_ns_solver.hpp b/include/unsteady_ns_solver.hpp index c579b57a..81c9636a 100644 --- a/include/unsteady_ns_solver.hpp +++ b/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 &basis_tags) override; diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 9b16800c..977e540e 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -25,15 +25,17 @@ AdvDiffSolver::AdvDiffSolver() load_flow = config.GetOption("adv-diff/load_flow", false); if (save_flow || load_flow) flow_file = config.GetRequiredOption("adv-diff/flow_file"); + + flow_solver_type = config.GetOption("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::functionGetNumVar() : 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; } \ No newline at end of file diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index b78b0f4c..fdbf57ee 100644 --- a/src/linelast_solver.cpp +++ b/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() diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index cec992a6..4a6b502b 100644 --- a/src/main_workflow.cpp +++ b/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(); diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index d946c710..915a2184 100644 --- a/src/multiblock_solver.cpp +++ b/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) diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index 948ed837..e0f514b9 100644 --- a/src/poisson_solver.cpp +++ b/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() diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index 3edc51f3..e86f580d 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -536,10 +536,10 @@ void MFEMROMHandler::Solve(BlockVector &rhs, BlockVector &sol) { assert(operator_loaded); - int maxIter = config.GetOption("solver/max_iter", 10000); - double rtol = config.GetOption("solver/relative_tolerance", 1.e-15); - double atol = config.GetOption("solver/absolute_tolerance", 1.e-15); - int print_level = config.GetOption("solver/print_level", 0); + int maxIter = config.GetOption("rom_solver/max_iter", 10000); + double rtol = config.GetOption("rom_solver/relative_tolerance", 1.e-15); + double atol = config.GetOption("rom_solver/absolute_tolerance", 1.e-15); + int print_level = config.GetOption("rom_solver/print_level", 0); std::string prec_str = config.GetOption("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("solver/use_restart", false); - double amp = config.GetOption("solver/initial_random_perturbation", 1.0e-5); + bool use_restart = config.GetOption("rom_solver/use_restart", false); + double amp = config.GetOption("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("solver/max_iter", 100); - double rtol = config.GetOption("solver/relative_tolerance", 1.e-10); - double atol = config.GetOption("solver/absolute_tolerance", 1.e-10); - int print_level = config.GetOption("solver/print_level", 0); + int maxIter = config.GetOption("rom_solver/max_iter", 100); + double rtol = config.GetOption("rom_solver/relative_tolerance", 1.e-10); + double atol = config.GetOption("rom_solver/absolute_tolerance", 1.e-10); + int print_level = config.GetOption("rom_solver/print_level", 0); - int jac_maxIter = config.GetOption("solver/jacobian/max_iter", 10000); - double jac_rtol = config.GetOption("solver/jacobian/relative_tolerance", 1.e-10); - double jac_atol = config.GetOption("solver/jacobian/absolute_tolerance", 1.e-10); - int jac_print_level = config.GetOption("solver/jacobian/print_level", -1); + int jac_maxIter = config.GetOption("rom_solver/jacobian/max_iter", 10000); + double jac_rtol = config.GetOption("rom_solver/jacobian/relative_tolerance", 1.e-10); + double jac_atol = config.GetOption("rom_solver/jacobian/absolute_tolerance", 1.e-10); + int jac_print_level = config.GetOption("rom_solver/jacobian/print_level", -1); std::string prec_str = config.GetOption("model_reduction/preconditioner", "none"); if (prec_str != "none") assert(prec); diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index 2ff664a1..4463faed 100644 --- a/src/stokes_solver.cpp +++ b/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; } -} \ No newline at end of file +} diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp index bda721fc..4d5dc739 100644 --- a/src/unsteady_ns_solver.cpp +++ b/src/unsteady_ns_solver.cpp @@ -290,7 +290,7 @@ void UnsteadyNSSolver::SaveVisualization(const int step, const double time) MultiBlockSolver::SaveVisualization(step, time); } -void UnsteadyNSSolver::SetParameterizedProblem(ParameterizedProblem *problem) +bool UnsteadyNSSolver::SetParameterizedProblem(ParameterizedProblem *problem) { SteadyNSSolver::SetParameterizedProblem(problem); @@ -313,7 +313,8 @@ void UnsteadyNSSolver::SetParameterizedProblem(ParameterizedProblem *problem) { u_ic = new VectorConstantCoefficient(zero_vel); p_ic = new VectorConstantCoefficient(zero_pres); - return; + // UnsteadyNSSolver does not fail in SetParameterizedProblem. + return true; } if (problem->ic_ptr[0]) @@ -325,6 +326,8 @@ void UnsteadyNSSolver::SetParameterizedProblem(ParameterizedProblem *problem) p_ic = new VectorFunctionCoefficient(1, problem->ic_ptr[1]); else p_ic = new VectorConstantCoefficient(zero_pres); + // UnsteadyNSSolver does not fail in SetParameterizedProblem. + return true; } BlockVector* UnsteadyNSSolver::PrepareSnapshots(std::vector &basis_tags) diff --git a/utils/gmsh2mfem.cpp b/utils/gmsh2mfem.cpp index 40122537..3eb5175c 100644 --- a/utils/gmsh2mfem.cpp +++ b/utils/gmsh2mfem.cpp @@ -13,6 +13,7 @@ using namespace mfem; // void trans(const Vector &x, Vector &p); static int order_ = 3; static bool force_nc_ = false; +static int ser_ref_levels = 0; int main(int argc, char *argv[]) { @@ -26,6 +27,8 @@ int main(int argc, char *argv[]) "Order (polynomial degree) of the mesh elements."); args.AddOption(&force_nc_, "-fnc", "--force-non-conforming", "-nfnc", "--noforce-non-conforming", "Sets whether to force the output mesh to be nonconforming. Default behavior is no enforcing."); + args.AddOption(&ser_ref_levels, "-rs", "--refine-serial", + "Number of times to refine the mesh uniformly in serial."); args.Parse(); if (!args.Good()) @@ -37,6 +40,12 @@ int main(int argc, char *argv[]) Mesh mesh(meshFileString); + // Refine the mesh if desired + for (int lev = 0; lev < ser_ref_levels; lev++) + { + mesh.UniformRefinement(); + } + // Promote to high order mesh if (order_ >1) mesh.SetCurvature(order_, true, 2, Ordering::byVDIM);