From 75ab3c0acc9295dd7bcdb80ee993605f62c25c9b Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 19 Dec 2025 10:51:13 -0800 Subject: [PATCH 01/27] advdiff: add option for steady-ns solver. --- include/advdiff_solver.hpp | 5 +++-- src/advdiff_solver.cpp | 46 +++++++++++++++++++++----------------- 2 files changed, 29 insertions(+), 22 deletions(-) diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index 62a39533..965c4d06 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 = ""; diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 9b16800c..2b8693fe 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() @@ -188,8 +190,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 +203,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]); @@ -240,38 +242,42 @@ void 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; + + 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_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(). */ } \ No newline at end of file From 502a8a2252438e5c2a1ee2f8a3eb38e98244427d Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Mon, 22 Dec 2025 09:36:22 -0800 Subject: [PATCH 02/27] error for unknown flow solver --- src/advdiff_solver.cpp | 2 ++ utils/gmsh2mfem.cpp | 9 +++++++++ 2 files changed, 11 insertions(+) diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 2b8693fe..890046e1 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -246,6 +246,8 @@ void AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem) 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(); 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); From 8b5aae2eef068b0833d070fdfa2213016c237f9d Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Mon, 22 Dec 2025 09:59:54 -0800 Subject: [PATCH 03/27] changing function signature of SetParameterizedProblem. --- include/advdiff_solver.hpp | 2 +- include/linelast_solver.hpp | 2 +- include/multiblock_solver.hpp | 2 +- include/poisson_solver.hpp | 2 +- include/stokes_solver.hpp | 2 +- include/unsteady_ns_solver.hpp | 2 +- src/advdiff_solver.cpp | 3 ++- src/linelast_solver.cpp | 3 ++- src/main_workflow.cpp | 7 ++++++- src/multiblock_solver.cpp | 3 ++- src/poisson_solver.cpp | 3 ++- src/stokes_solver.cpp | 2 +- src/unsteady_ns_solver.cpp | 5 +++-- 13 files changed, 24 insertions(+), 14 deletions(-) diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index 965c4d06..2bcb45cd 100644 --- a/include/advdiff_solver.hpp +++ b/include/advdiff_solver.hpp @@ -53,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; 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 890046e1..aa09bb42 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -171,12 +171,13 @@ void AdvDiffSolver::SetFlowAtSubdomain(std::functiongeneral_vector_ptr[0])); } + return true; } void LinElastSolver::ProjectOperatorOnReducedBasis() diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index cec992a6..3e515236 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -165,7 +165,12 @@ void GenerateSamples(MPI_Comm comm) test->InitROMHandler(); problem->SetSingleRun(); - test->SetParameterizedProblem(problem); + bool param_set = test->SetParameterizedProblem(problem); + if (!param_set) + { + sample_generator->SetSampleParams(s); + param_set = test->SetParameterizedProblem(problem); + } int file_idx = s + sample_generator->GetFileOffset(); const std::string visual_path = sample_generator->GetSamplePath(file_idx, test->GetVisualizationPrefix()); diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index d946c710..00161818 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,7 @@ void MultiBlockSolver::SetParameterizedProblem(ParameterizedProblem *problem) bdr_type[b] = problem->bdr_type[idx]; } + return true; } void MultiBlockSolver::SaveSolution(std::string filename) diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index 948ed837..9c3fbada 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,7 @@ void PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem) AddRHSFunction(*(problem->scalar_rhs_ptr)); else AddRHSFunction(0.0); + return true; } void PoissonSolver::SetMUMPSSolver() diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index 2ff664a1..69457f7d 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); diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp index bda721fc..24d64119 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,7 @@ void UnsteadyNSSolver::SetParameterizedProblem(ParameterizedProblem *problem) { u_ic = new VectorConstantCoefficient(zero_vel); p_ic = new VectorConstantCoefficient(zero_pres); - return; + return true; } if (problem->ic_ptr[0]) @@ -325,6 +325,7 @@ void UnsteadyNSSolver::SetParameterizedProblem(ParameterizedProblem *problem) p_ic = new VectorFunctionCoefficient(1, problem->ic_ptr[1]); else p_ic = new VectorConstantCoefficient(zero_pres); + return true; } BlockVector* UnsteadyNSSolver::PrepareSnapshots(std::vector &basis_tags) From cb0db4d00acf6a56b8838356ca6a0611d940c826 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Mon, 22 Dec 2025 10:11:21 -0800 Subject: [PATCH 04/27] resample if SetParameterizedProblem fails. --- include/advdiff_solver.hpp | 2 +- src/advdiff_solver.cpp | 13 ++++++++----- src/linelast_solver.cpp | 1 + src/main_workflow.cpp | 6 +++--- src/multiblock_solver.cpp | 1 + src/poisson_solver.cpp | 1 + src/unsteady_ns_solver.cpp | 2 ++ 7 files changed, 17 insertions(+), 9 deletions(-) diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index 2bcb45cd..c39691d7 100644 --- a/include/advdiff_solver.hpp +++ b/include/advdiff_solver.hpp @@ -61,7 +61,7 @@ friend class ParameterizedProblem; void SetMUMPSSolver() override; private: - void GetFlowField(ParameterizedProblem *flow_problem); + bool GetFlowField(ParameterizedProblem *flow_problem); }; #endif diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index aa09bb42..977e540e 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -173,11 +173,12 @@ void AdvDiffSolver::SetFlowAtSubdomain(std::functionSetOperator(*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"); @@ -267,7 +268,7 @@ void AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem) flow_solver->BuildOperators(); flow_solver->SetupBCOperators(); flow_solver->Assemble(); - flow_solver->Solve(); + flow_loaded = flow_solver->Solve(); } if (save_flow && (!flow_loaded)) @@ -283,4 +284,6 @@ void AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem) and it requires the grid function for its lifetime. 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 6042d893..fdbf57ee 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -483,6 +483,7 @@ bool LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem) { SetupIC(*(problem->general_vector_ptr[0])); } + // LinElastSolver does not fail in SetParameterizedProblem. return true; } diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index 3e515236..f7972efe 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -432,7 +432,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(); @@ -509,7 +509,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(); @@ -726,7 +726,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 00161818..915a2184 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -502,6 +502,7 @@ bool MultiBlockSolver::SetParameterizedProblem(ParameterizedProblem *problem) bdr_type[b] = problem->bdr_type[idx]; } + // MultiBlockSolver does not fail in SetParameterizedProblem. return true; } diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index 9c3fbada..e0f514b9 100644 --- a/src/poisson_solver.cpp +++ b/src/poisson_solver.cpp @@ -582,6 +582,7 @@ bool PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem) AddRHSFunction(*(problem->scalar_rhs_ptr)); else AddRHSFunction(0.0); + // PoissonSolver does not fail in SetParameterizedProblem. return true; } diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp index 24d64119..4d5dc739 100644 --- a/src/unsteady_ns_solver.cpp +++ b/src/unsteady_ns_solver.cpp @@ -313,6 +313,7 @@ bool UnsteadyNSSolver::SetParameterizedProblem(ParameterizedProblem *problem) { u_ic = new VectorConstantCoefficient(zero_vel); p_ic = new VectorConstantCoefficient(zero_pres); + // UnsteadyNSSolver does not fail in SetParameterizedProblem. return true; } @@ -325,6 +326,7 @@ bool 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; } From bdbccbe41120f603485e26230f2f69936ea18645 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Mon, 22 Dec 2025 10:17:31 -0800 Subject: [PATCH 05/27] error message if setparameterizedproblem fails for non-random sampling. --- src/main_workflow.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index f7972efe..e3d3d414 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -168,6 +168,8 @@ void GenerateSamples(MPI_Comm comm) 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"); sample_generator->SetSampleParams(s); param_set = test->SetParameterizedProblem(problem); } From 6d045e7181eec888dd35af12cf4ecf121b3cf030 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Mon, 22 Dec 2025 10:22:19 -0800 Subject: [PATCH 06/27] fix stokes solver. --- src/stokes_solver.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index 69457f7d..4463faed 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -1143,6 +1143,9 @@ bool 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 +} From 8397e4d18b41914ff74e15f9a765b688e3cd0ee0 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Fri, 26 Dec 2025 23:02:11 -0800 Subject: [PATCH 07/27] fix resampling --- src/main_workflow.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index e3d3d414..4a6b502b 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -170,7 +170,15 @@ void GenerateSamples(MPI_Comm comm) { 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); } From 95154aaf215c5c465658527934b692300b083858 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Fri, 26 Dec 2025 23:02:39 -0800 Subject: [PATCH 08/27] rom solver follows different options than fom solver. --- src/rom_handler.cpp | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) 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); From aaacd90fafcf66099093bd9b3a894197e7dba1e4 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Mon, 12 Jan 2026 13:55:45 -0800 Subject: [PATCH 09/27] use different order for flow solver --- include/advdiff_solver.hpp | 1 + src/advdiff_solver.cpp | 11 ++++++++++- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index c39691d7..c761d924 100644 --- a/include/advdiff_solver.hpp +++ b/include/advdiff_solver.hpp @@ -27,6 +27,7 @@ friend class ParameterizedProblem; flow solver to obtain the prescribed velocity field. both StokesSolver / SteadyNSSolver can be used. */ std::string flow_solver_type = ""; + int flow_solver_order = -1; StokesSolver *flow_solver = NULL; bool load_flow = false; bool save_flow = false; diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 977e540e..c28ffed7 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -26,7 +26,8 @@ AdvDiffSolver::AdvDiffSolver() if (save_flow || load_flow) flow_file = config.GetRequiredOption("adv-diff/flow_file"); - flow_solver_type = config.GetOption("adv-diff/flow_solver_type", "stokes"); + flow_solver_type = config.GetOption("adv-diff/flow_solver/type", "stokes"); + flow_solver_order = config.GetOption("adv-diff/flow_solver/order", order); } AdvDiffSolver::~AdvDiffSolver() @@ -244,6 +245,10 @@ 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"); + // Temporarily change order option just for flow solver setup. + if (flow_solver_order != order) + config.SetOption("discretization/order", flow_solver_order); + if (flow_solver_type == "stokes") flow_solver = new StokesSolver; else if (flow_solver_type == "steady-ns") @@ -251,6 +256,10 @@ bool AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem) else mfem_error("AdvDiffSolver::GetFlowField - Unknown flow solver type!\n"); + // Bring back the original value for order option. + if (flow_solver_order != order) + config.SetOption("discretization/order", order); + flow_solver->InitVariables(); if (use_rom) flow_solver->InitROMHandler(); flow_solver->SetSolutionSaveMode(save_flow); From f19424871dd2f297ec2d07af8c85e7f71be1751d Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 15 Jan 2026 11:45:54 -0800 Subject: [PATCH 10/27] ComponentTopologyHandler subset constructor. --- include/component_topology_handler.hpp | 15 ++ src/component_topology_handler.cpp | 261 +++++++++++++++++++++++++ 2 files changed, 276 insertions(+) diff --git a/include/component_topology_handler.hpp b/include/component_topology_handler.hpp index c90222ec..294706c5 100644 --- a/include/component_topology_handler.hpp +++ b/include/component_topology_handler.hpp @@ -83,6 +83,13 @@ class ComponentTopologyHandler : public TopologyHandler std::unordered_map vtx2to1; // vertex mapping from component 2 to component 1. Array2D 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,14 @@ 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. + // Boundary attributes of subset meshes facing outside are set as global boundary. + ComponentTopologyHandler(ComponentTopologyHandler* global, const int i0, const int j0, const int N, const int M); + virtual ~ComponentTopologyHandler(); // access diff --git a/src/component_topology_handler.cpp b/src/component_topology_handler.cpp index 25a2a863..d47c5c5b 100644 --- a/src/component_topology_handler.cpp +++ b/src/component_topology_handler.cpp @@ -1008,4 +1008,265 @@ bool ComponentTopologyHandler::ComponentBdrAttrCheck(Mesh *comp) } return success; +} + +ComponentTopologyHandler::ComponentTopologyHandler( + ComponentTopologyHandler* global, const int i0, const int j0, + const int N, const int M) + : TopologyHandler(TopologyHandlerMode::COMPONENT), + verbose(global->verbose), write_ports(false), vtx_gap_thrs(global->vtx_gap_thrs), + tf_ptr(global->tf_ptr), inv_tf_ptr(global->inv_tf_ptr) +{ + // assert N from the N x N array + assert(N * N == numSub); + + // Validate input parameters + assert(i0 >= 0 && j0 >= 0); + assert(i0 + M <= N && j0 + M <= N); + + // Set up mesh arrays for the M x M subset + int new_numSub = M * M; + numSub = new_numSub; + + // Copy basic topology information + dim = global->dim; + dd_mode = global->dd_mode; + num_comp = global->num_comp; + comp_names = global->comp_names; + + // Copy components (reference meshes) + components.SetSize(num_comp); + for (int c = 0; c < num_comp; c++) + components[c] = new BlockMesh(*global->components[c]); + + sub_composition.SetSize(new_numSub); + sub_composition = 0; + mesh_types.SetSize(new_numSub); + mesh_comp_idx.SetSize(new_numSub); + mesh_configs.SetSize(new_numSub); + meshes.SetSize(new_numSub); + bdr_c2g.SetSize(new_numSub); + + Array comp_included(global->num_comp); + comp_included = 0; + + // Map from subset mesh index to original mesh index + Array orig_to_subset(global->numSub); + orig_to_subset = -1; + Array subset_to_orig(new_numSub); + subset_to_orig = -1; + + // Copy meshes and configurations for the subset + for (int mi = 0; mi < M; mi++) + { + for (int mj = 0; mj < M; mj++) + { + int subset_idx = mi * M + mj; + int orig_idx = (i0 + mi) * N + (j0 + mj); + + subset_to_orig[subset_idx] = orig_idx; + orig_to_subset[orig_idx] = subset_idx; + + // Track which components are actually used + int orig_comp_type = global->mesh_types[orig_idx]; + comp_included[orig_comp_type] = 1; + + // Copy mesh type (will be remapped later) + mesh_types[subset_idx] = orig_comp_type; + + // Copy mesh configuration + for (int d = 0; d < 3; d++) + { + mesh_configs[subset_idx].trans[d] = global->mesh_configs[orig_idx].trans[d]; + mesh_configs[subset_idx].rotate[d] = global->mesh_configs[orig_idx].rotate[d]; + } + + // Copy the mesh (as-is, no transformation needed since it's already transformed) + meshes[subset_idx] = new Mesh(*global->meshes[orig_idx]); + + // Copy boundary attribute mapping for this mesh + bdr_c2g[subset_idx] = new Array(*global->bdr_c2g[orig_idx]); + } + } + + // Now update num_comp, components, comp_names based on actually used components + Array old_to_new_comp(global->num_comp); + old_to_new_comp = -1; + num_comp = 0; + + for (int c = 0; c < global->num_comp; c++) + { + if (comp_included[c]) + { + old_to_new_comp[c] = num_comp; + num_comp++; + } + } + + // Resize and populate the component arrays + components.SetSize(num_comp); + comp_names.resize(num_comp); + comp_name2idx.clear(); + + for (int c = 0; c < global->num_comp; c++) + { + if (comp_included[c]) + { + int new_c = old_to_new_comp[c]; + components[new_c] = new BlockMesh(*global->components[c]); + comp_names[new_c] = global->comp_names[c]; + comp_name2idx[comp_names[new_c]] = new_c; + } + } + + // Remap mesh_types to new component indices + for (int m = 0; m < new_numSub; m++) + { + mesh_types[m] = old_to_new_comp[mesh_types[m]]; + } + + // Calculate sub_composition and mesh_comp_idx based on new component indices + sub_composition = 0; + for (int m = 0; m < new_numSub; m++) + { + int comp_type = mesh_types[m]; + mesh_comp_idx[m] = sub_composition[comp_type]; + sub_composition[comp_type] += 1; + } + + // Determine which ports are internal (both meshes in subset) + // and which are external (at least one mesh outside subset) + Array port_included(global->num_ports); + port_included = 0; + + Array new_port_infos(0); + Array new_port_types(0); + Array ref_port_included(global->num_ref_ports); + ref_port_included = 0; + + // Process each port to determine if it's internal or external + for (int p = 0; p < global->num_ports; p++) + { + int mesh1 = global->port_infos[p].Mesh1; + int mesh2 = global->port_infos[p].Mesh2; + int attr1 = global->port_infos[p].Attr1; + int attr2 = global->port_infos[p].Attr2; + int port_attr = global->port_infos[p].PortAttr; + int port_type = global->port_types[p]; + + // Check if both meshes are in the subset + bool mesh1_in_subset = (orig_to_subset[mesh1] >= 0); + bool mesh2_in_subset = (orig_to_subset[mesh2] >= 0); + + if (mesh1_in_subset && mesh2_in_subset) + { + // Internal port - include it + port_included[p] = 1; + ref_port_included[port_type] = 1; + + int new_mesh1 = orig_to_subset[mesh1]; + int new_mesh2 = orig_to_subset[mesh2]; + + // Add to new port infos with updated mesh indices + PortInfo new_info = global->port_infos[p]; + new_info.Mesh1 = new_mesh1; + new_info.Mesh2 = new_mesh2; + new_port_infos.Append(new_info); + new_port_types.Append(port_type); // Still using old port_type index, will remap later + } + else + { + // External port - exclude it + // The boundary attributes of meshes facing outside become global boundary + port_included[p] = 0; + } + } + + // Create mapping from old to new ref_port indices + Array old_to_new_refport(global->num_ref_ports); + old_to_new_refport = -1; + num_ref_ports = 0; + + for (int rp = 0; rp < global->num_ref_ports; rp++) + { + if (ref_port_included[rp]) + { + old_to_new_refport[rp] = num_ref_ports; + num_ref_ports++; + } + } + + // Resize and populate the reference port arrays + ref_ports.SetSize(num_ref_ports); + port_names.resize(num_ref_ports); + port_dicts.SetSize(num_ref_ports); + ref_interfaces.SetSize(num_ref_ports); + port_name2idx.clear(); + + for (int rp = 0; rp < global->num_ref_ports; rp++) + { + if (ref_port_included[rp]) + { + int new_rp = old_to_new_refport[rp]; + + // Copy PortData and update component indices + ref_ports[new_rp] = new PortData(*global->ref_ports[rp]); + ref_ports[new_rp]->Component1 = old_to_new_comp[global->ref_ports[rp]->Component1]; + ref_ports[new_rp]->Component2 = old_to_new_comp[global->ref_ports[rp]->Component2]; + + // Copy port name + port_names[new_rp] = global->port_names[rp]; + port_name2idx[port_names[new_rp]] = new_rp; + + // Copy port dictionary + port_dicts[new_rp] = new YAML::Node(*global->port_dicts[rp]); + + // Copy reference interface + ref_interfaces[new_rp] = new Array(*global->ref_interfaces[rp]); + } + } + + // Remap port_types in new_port_types to new reference port indices + for (int p = 0; p < new_port_types.Size(); p++) + { + new_port_types[p] = old_to_new_refport[new_port_types[p]]; + } + + // Assign the updated arrays + port_infos = new_port_infos; + port_types = new_port_types; + num_ports = port_infos.Size(); + + // Setup ports, as all reference port infos are updated. + SetupPorts(); + + // Update bdr_attributes based on actual boundary attributes in the subset + // Start by collecting all boundary attributes from bdr_c2g + std::set bdr_attr_set; + + for (int m = 0; m < new_numSub; m++) + { + for (int i = 0; i < bdr_c2g[m]->Size(); i++) + { + int attr = (*bdr_c2g[m])[i]; + bdr_attr_set.insert(attr); + } + } + + // Exclude port attributes (these are internal interfaces, not true boundaries) + for (int p = 0; p < port_infos.Size(); p++) + { + bdr_attr_set.erase(port_infos[p].PortAttr); + } + + // Convert set to array + bdr_attributes.SetSize(bdr_attr_set.size()); + int idx = 0; + for (std::set::iterator it = bdr_attr_set.begin(); it != bdr_attr_set.end(); ++it) + { + bdr_attributes[idx] = *it; + idx++; + } + + bdr_attributes.Sort(); } \ No newline at end of file From 1ed5556a3ab43b0b6031a74198a4604dfc94deb3 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 15 Jan 2026 11:56:43 -0800 Subject: [PATCH 11/27] MultiBlockSolver can take a-priori defined TopologyHandler. --- include/advdiff_solver.hpp | 2 +- include/linelast_solver.hpp | 2 +- include/multiblock_solver.hpp | 2 +- include/poisson_solver.hpp | 2 +- include/steady_ns_solver.hpp | 2 +- include/stokes_solver.hpp | 2 +- include/unsteady_ns_solver.hpp | 2 +- src/advdiff_solver.cpp | 4 ++-- src/linelast_solver.cpp | 4 ++-- src/multiblock_solver.cpp | 36 ++++++++++++++++++---------------- src/poisson_solver.cpp | 4 ++-- src/steady_ns_solver.cpp | 4 ++-- src/stokes_solver.cpp | 4 ++-- src/unsteady_ns_solver.cpp | 4 ++-- 14 files changed, 38 insertions(+), 36 deletions(-) diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index c761d924..9990c681 100644 --- a/include/advdiff_solver.hpp +++ b/include/advdiff_solver.hpp @@ -41,7 +41,7 @@ friend class ParameterizedProblem; Array global_flow_visual; public: - AdvDiffSolver(); + AdvDiffSolver(TopologyHandler *input_topol_handler=NULL); virtual ~AdvDiffSolver(); diff --git a/include/linelast_solver.hpp b/include/linelast_solver.hpp index e499aae4..c3207010 100644 --- a/include/linelast_solver.hpp +++ b/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(); diff --git a/include/multiblock_solver.hpp b/include/multiblock_solver.hpp index d7f3ca4c..8ae055c4 100644 --- a/include/multiblock_solver.hpp +++ b/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(); diff --git a/include/poisson_solver.hpp b/include/poisson_solver.hpp index 115dd1ee..32216c4d 100644 --- a/include/poisson_solver.hpp +++ b/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(); diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index c5bb02c8..5c35d021 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -166,7 +166,7 @@ friend class SteadyNSOperator; NewtonSolver *newton_solver = NULL; public: - SteadyNSSolver(); + SteadyNSSolver(TopologyHandler *input_topol_handler=NULL); virtual ~SteadyNSSolver(); diff --git a/include/stokes_solver.hpp b/include/stokes_solver.hpp index 614b41b9..d5e16b65 100644 --- a/include/stokes_solver.hpp +++ b/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(); diff --git a/include/unsteady_ns_solver.hpp b/include/unsteady_ns_solver.hpp index 81c9636a..eb60539c 100644 --- a/include/unsteady_ns_solver.hpp +++ b/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(); diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index c28ffed7..89670b5e 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -10,8 +10,8 @@ using namespace std; using namespace mfem; -AdvDiffSolver::AdvDiffSolver() - : PoissonSolver(), flow_visual(0), flow_fes(0), global_flow_visual(0) +AdvDiffSolver::AdvDiffSolver(TopologyHandler *input_topol_handler) + : PoissonSolver(input_topol_handler), flow_visual(0), flow_fes(0), global_flow_visual(0) { // ConvectionIntegrator does not support L2 space. assert(!full_dg); diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index fdbf57ee..696a847b 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -10,8 +10,8 @@ using namespace std; using namespace mfem; -LinElastSolver::LinElastSolver() - : MultiBlockSolver() +LinElastSolver::LinElastSolver(TopologyHandler *input_topol_handler) + : MultiBlockSolver(input_topol_handler) { alpha = config.GetOption("discretization/interface/alpha", -1.0); kappa = config.GetOption("discretization/interface/kappa", (order + 1) * (order + 1)); diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index 915a2184..ef28c1a3 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -10,7 +10,8 @@ using namespace std; using namespace mfem; -MultiBlockSolver::MultiBlockSolver() +MultiBlockSolver::MultiBlockSolver(TopologyHandler *input_topol_handler) + : topol_handler(input_topol_handler) { /* TODO(kevin): this is a boilerplate for parallel POD/EQP training. @@ -22,24 +23,25 @@ MultiBlockSolver::MultiBlockSolver() ParseInputs(); TopologyData topol_data; - switch (topol_mode) - { - case TopologyHandlerMode::SUBMESH: - { - topol_handler = new SubMeshTopologyHandler(); - break; - } - case TopologyHandlerMode::COMPONENT: + if (topol_handler == NULL) + switch (topol_mode) { - topol_handler = new ComponentTopologyHandler(); - break; - } - default: - { - mfem_error("Unknown topology handler mode!\n"); - break; + case TopologyHandlerMode::SUBMESH: + { + topol_handler = new SubMeshTopologyHandler(); + break; + } + case TopologyHandlerMode::COMPONENT: + { + topol_handler = new ComponentTopologyHandler(); + break; + } + default: + { + mfem_error("Unknown topology handler mode!\n"); + break; + } } - } topol_handler->ExportInfo(meshes, topol_data); // Receive topology info diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index e0f514b9..bcaccb8b 100644 --- a/src/poisson_solver.cpp +++ b/src/poisson_solver.cpp @@ -10,8 +10,8 @@ using namespace std; using namespace mfem; -PoissonSolver::PoissonSolver() - : MultiBlockSolver() +PoissonSolver::PoissonSolver(TopologyHandler *input_topol_handler) + : MultiBlockSolver(input_topol_handler) { sigma = config.GetOption("discretization/interface/sigma", -1.0); kappa = config.GetOption("discretization/interface/kappa", (order + 1) * (order + 1)); diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index 7a3ae61a..560bf8bd 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -416,8 +416,8 @@ void SteadyNSEQPROM::AddVel(const Vector &y_u, Vector &y) const SteadyNSSolver */ -SteadyNSSolver::SteadyNSSolver() - : StokesSolver() +SteadyNSSolver::SteadyNSSolver(TopologyHandler *input_topol_handler) + : StokesSolver(input_topol_handler) { nonlinear_mode = true; diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index 4463faed..f663e729 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -13,8 +13,8 @@ using namespace std; using namespace mfem; -StokesSolver::StokesSolver() - : MultiBlockSolver(), minus_one(-1.0) +StokesSolver::StokesSolver(TopologyHandler *input_topol_handler) + : MultiBlockSolver(input_topol_handler), minus_one(-1.0) { nu = config.GetOption("stokes/nu", 1.0); nu_coeff = new ConstantCoefficient(nu); diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp index 4d5dc739..d70fa040 100644 --- a/src/unsteady_ns_solver.cpp +++ b/src/unsteady_ns_solver.cpp @@ -12,8 +12,8 @@ using namespace mfem; UnsteadyNSSolver */ -UnsteadyNSSolver::UnsteadyNSSolver() - : SteadyNSSolver(), timer("UnsteadyNSSolver") +UnsteadyNSSolver::UnsteadyNSSolver(TopologyHandler *input_topol_handler) + : SteadyNSSolver(input_topol_handler), timer("UnsteadyNSSolver") { nonlinear_mode = true; From 0ebb1eb83b2f004ef76f82209edebaadcb00fb6d Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 15 Jan 2026 15:48:17 -0800 Subject: [PATCH 12/27] fix linelast BC setup to be consistent with other solver. option fix for steady ns solver test. --- include/linelast_solver.hpp | 3 ++ src/linelast_solver.cpp | 69 ++++++++++++++++++++++++----- test/inputs/steady_ns.base.yml | 7 +++ test/inputs/steady_ns.component.yml | 7 +++ 4 files changed, 76 insertions(+), 10 deletions(-) diff --git a/include/linelast_solver.hpp b/include/linelast_solver.hpp index c3207010..6e4d88af 100644 --- a/include/linelast_solver.hpp +++ b/include/linelast_solver.hpp @@ -53,6 +53,8 @@ class LinElastSolver : public MultiBlockSolver double lambda = 1.0; double mu = 1.0; + Vector zero; + // Initial positions VectorCoefficient *init_x = NULL; @@ -86,6 +88,7 @@ class LinElastSolver : public MultiBlockSolver virtual void SetupBCVariables() override; virtual void SetupIC(std::function F); virtual void AddBCFunction(std::function F, const int battr = -1); + virtual void AddBCFunction(const Vector &F, const int battr = -1); virtual void AddRHSFunction(std::function F); virtual bool BCExistsOnBdr(const int &global_battr_idx); virtual void SetupBCOperators(); diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index 696a847b..d2a8d83b 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -40,6 +40,10 @@ LinElastSolver::LinElastSolver(TopologyHandler *input_topol_handler) { fes[m] = new FiniteElementSpace(meshes[m], fec[0], udim); } + + // For homogeneous Dirichlet condition + zero.SetSize(dim); + zero = 0.0; } LinElastSolver::~LinElastSolver() @@ -178,6 +182,7 @@ void LinElastSolver::SetupRHSBCOperators() switch (bdr_type[b]) { + case BoundaryType::ZERO: case BoundaryType::DIRICHLET: bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator( *bdr_coeffs[b], *lambda_c[m], *mu_c[m], alpha, kappa), *bdr_markers[b]); @@ -406,6 +411,26 @@ void LinElastSolver::AddBCFunction(std::function 0); + + if (battr > 0) + { + int idx = global_bdr_attributes.Find(battr); + if (idx < 0) + { + std::string msg = "battr " + std::to_string(battr) + " is not in global boundary attributes. skipping this boundary condition.\n"; + mfem_warning(msg.c_str()); + return; + } + bdr_coeffs[idx] = new VectorConstantCoefficient(F); + } + else + for (int k = 0; k < bdr_coeffs.Size(); k++) + bdr_coeffs[k] = new VectorConstantCoefficient(F); +} + void LinElastSolver::AddRHSFunction(std::function F) { rhs_coeffs.Append(new VectorFunctionCoefficient(dim, F)); @@ -431,10 +456,13 @@ void LinElastSolver::SetupDomainBCOperators() continue; if (!BCExistsOnBdr(b)) continue; + /* Neumann condition is implicitly enforced by weak form */ + if (bdr_type[b] == BoundaryType::NEUMANN) + continue; - if (bdr_type[b] == BoundaryType::DIRICHLET) - as[m]->AddBdrFaceIntegrator(new DGElasticityIntegrator( - *(lambda_c[m]), *(mu_c[m]), alpha, kappa), *(bdr_markers[b])); + /* for DIRICHLET and ZERO */ + as[m]->AddBdrFaceIntegrator(new DGElasticityIntegrator( + *(lambda_c[m]), *(mu_c[m]), alpha, kappa), *(bdr_markers[b])); } } } @@ -464,13 +492,34 @@ bool LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem) // Set BCs, the switch on BC type is done inside SetupRHSBCOperators for (int b = 0; b < problem->battr.Size(); b++) { - /* Dirichlet bc requires a function specified, even for zero. */ - if (problem->bdr_type[b] == BoundaryType::DIRICHLET) - assert(problem->vector_bdr_ptr[b]); - - /* Neumann bc does not require a function specified for zero */ - if (problem->vector_bdr_ptr[b]) - AddBCFunction(*(problem->vector_bdr_ptr[b]), problem->battr[b]); + switch (problem->bdr_type[b]) + { + case BoundaryType::DIRICHLET: + { + assert(problem->vector_bdr_ptr[b]); + AddBCFunction(*(problem->vector_bdr_ptr[b]), problem->battr[b]); + break; + } + case BoundaryType::NEUMANN: + { + if (problem->vector_bdr_ptr[b]) + AddBCFunction(*(problem->vector_bdr_ptr[b]), problem->battr[b]); + break; + } + default: + case BoundaryType::ZERO: + { AddBCFunction(zero, problem->battr[b]); break; } + } + } + + /* + for all undefined boundary attributes, + set homogeneous Neumann condition. + */ + for (int k = 0; k < bdr_type.Size(); k++) + { + if (bdr_type[k] == BoundaryType::NUM_BDR_TYPE) + bdr_type[k] = BoundaryType::NEUMANN; } // Set RHS diff --git a/test/inputs/steady_ns.base.yml b/test/inputs/steady_ns.base.yml index 6265572d..bed575d9 100644 --- a/test/inputs/steady_ns.base.yml +++ b/test/inputs/steady_ns.base.yml @@ -23,6 +23,13 @@ solver: relative_tolerance: 1.0e-15 absolute_tolerance: 1.0e-15 +rom_solver: + relative_tolerance: 1.0e-15 + absolute_tolerance: 1.0e-15 + jacobian: + relative_tolerance: 1.0e-15 + absolute_tolerance: 1.0e-15 + visualization: enabled: false unified_paraview: true diff --git a/test/inputs/steady_ns.component.yml b/test/inputs/steady_ns.component.yml index d2f6b174..3ed86d96 100644 --- a/test/inputs/steady_ns.component.yml +++ b/test/inputs/steady_ns.component.yml @@ -47,6 +47,13 @@ solver: relative_tolerance: 1.0e-15 absolute_tolerance: 1.0e-15 +rom_solver: + relative_tolerance: 1.0e-15 + absolute_tolerance: 1.0e-15 + jacobian: + relative_tolerance: 1.0e-15 + absolute_tolerance: 1.0e-15 + visualization: enabled: false unified_paraview: true From a7568b9e1cbdb059b866482bc624fef2b26fe13f Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 15 Jan 2026 16:11:21 -0800 Subject: [PATCH 13/27] sanity check on bdr type definition. --- include/multiblock_solver.hpp | 9 +++++++++ src/linelast_solver.cpp | 1 + src/main_workflow.cpp | 2 ++ src/multiblock_solver.cpp | 7 ++++++- src/poisson_solver.cpp | 2 ++ src/stokes_solver.cpp | 2 ++ 6 files changed, 22 insertions(+), 1 deletion(-) diff --git a/include/multiblock_solver.hpp b/include/multiblock_solver.hpp index 8ae055c4..c23f88f3 100644 --- a/include/multiblock_solver.hpp +++ b/include/multiblock_solver.hpp @@ -262,6 +262,15 @@ friend class ParameterizedProblem; 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; + } + protected: virtual void AssembleROMMat(BlockMatrix &romMat); }; diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index d2a8d83b..acde4c92 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -438,6 +438,7 @@ void LinElastSolver::AddRHSFunction(std::functionSetSingleRun(); assert(test->SetParameterizedProblem(problem)); + assert(test->IsBdrTypeDefined()); + // TODO: there are skippable operations depending on rom/fom mode. test->BuildRHSOperators(); test->SetupRHSBCOperators(); diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index ef28c1a3..88ecab56 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -309,7 +309,12 @@ void MultiBlockSolver::AssembleROMMat(BlockMatrix &romMat) { int global_idx = global_bdr_attributes.Find((*bdr_c2g)[b]); if (global_idx < 0) continue; - if (!BCExistsOnBdr(global_idx)) continue; + /* + we do not depend on the existence of bc coefficients. + this is to support alternating Schwarz solver, + where internal boundaries would not have global bc coefficients. + */ + // if (!BCExistsOnBdr(global_idx)) continue; /* we assume only Neumann condition would not add an operator. */ if (bdr_type[global_idx] == BoundaryType::NEUMANN) diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index bcaccb8b..cdc50b7c 100644 --- a/src/poisson_solver.cpp +++ b/src/poisson_solver.cpp @@ -204,6 +204,8 @@ bool PoissonSolver::BCExistsOnBdr(const int &global_battr_idx) void PoissonSolver::SetupBCOperators() { + assert(IsBdrTypeDefined()); + SetupRHSBCOperators(); SetupDomainBCOperators(); diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index f663e729..3be19f43 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -303,6 +303,8 @@ bool StokesSolver::BCExistsOnBdr(const int &global_battr_idx) void StokesSolver::SetupBCOperators() { + assert(IsBdrTypeDefined()); + SetupRHSBCOperators(); SetupDomainBCOperators(); From 3e63c573e5beaff921ca91e2413f782724bec6f3 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 16 Jan 2026 09:23:19 -0800 Subject: [PATCH 14/27] subset constructor saves subset-to-original mesh mapping. --- include/component_topology_handler.hpp | 4 +++- src/component_topology_handler.cpp | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/include/component_topology_handler.hpp b/include/component_topology_handler.hpp index 294706c5..a2848182 100644 --- a/include/component_topology_handler.hpp +++ b/include/component_topology_handler.hpp @@ -146,7 +146,9 @@ class ComponentTopologyHandler : public TopologyHandler // Internal ports within the subset are included. // External ports connecting to meshes outside the subset are excluded. // Boundary attributes of subset meshes facing outside are set as global boundary. - ComponentTopologyHandler(ComponentTopologyHandler* global, const int i0, const int j0, const int N, const int M); + // Stores subset mesh to original mesh mapping. + ComponentTopologyHandler(ComponentTopologyHandler* global, const int i0, const int j0, + const int N, const int M, Array &subset_to_orig); virtual ~ComponentTopologyHandler(); diff --git a/src/component_topology_handler.cpp b/src/component_topology_handler.cpp index d47c5c5b..2ffd2e89 100644 --- a/src/component_topology_handler.cpp +++ b/src/component_topology_handler.cpp @@ -1012,7 +1012,7 @@ bool ComponentTopologyHandler::ComponentBdrAttrCheck(Mesh *comp) ComponentTopologyHandler::ComponentTopologyHandler( ComponentTopologyHandler* global, const int i0, const int j0, - const int N, const int M) + const int N, const int M, Array &subset_to_orig) : TopologyHandler(TopologyHandlerMode::COMPONENT), verbose(global->verbose), write_ports(false), vtx_gap_thrs(global->vtx_gap_thrs), tf_ptr(global->tf_ptr), inv_tf_ptr(global->inv_tf_ptr) @@ -1053,7 +1053,7 @@ ComponentTopologyHandler::ComponentTopologyHandler( // Map from subset mesh index to original mesh index Array orig_to_subset(global->numSub); orig_to_subset = -1; - Array subset_to_orig(new_numSub); + subset_to_orig.SetSize(new_numSub); subset_to_orig = -1; // Copy meshes and configurations for the subset From a6a9c566db30f7a2305958e50672af7fcae2e271 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 16 Jan 2026 09:24:24 -0800 Subject: [PATCH 15/27] StokesSolver::SetSubsetComplementaryFlux for subset incompressibility. --- include/stokes_solver.hpp | 7 ++ src/stokes_solver.cpp | 188 +++++++++++++++++++++++++++++++++++++- 2 files changed, 194 insertions(+), 1 deletion(-) diff --git a/include/stokes_solver.hpp b/include/stokes_solver.hpp index d5e16b65..54743794 100644 --- a/include/stokes_solver.hpp +++ b/include/stokes_solver.hpp @@ -216,6 +216,13 @@ friend class ParameterizedProblem; double ComputeBEIntegral(const FiniteElement &el, ElementTransformation &Tr, Coefficient &Q); void ComputeBEIntegral(const FiniteElement &el, ElementTransformation &Tr, VectorCoefficient &Q, Vector &result); + + // Assumes meshes are arranged in an N x N square array. + // Extracts an M x M subset starting at (i0, j0). + // to ensure incompressibility for the subset domain with all velocity dirichlet bc. + void SetSubsetComplementaryFlux(const int N, const int M, const int i0, const int j0, + const Array &subset_bdr_attributes, + const Array &subset_bdrtype); }; #endif diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index 3be19f43..67c19492 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -13,6 +13,28 @@ using namespace std; using namespace mfem; +namespace subset_flux +{ + double del_u; + Vector x0; + + void dir(const Vector &x, Vector &y) + { + const int dim = x.Size(); + y.SetSize(dim); + y = x; + + assert(x0.Size() == dim); + y -= x0; + } + + void flux(const Vector &x, Vector &y) + { + dir(x, y); + y *= del_u; + } +} + StokesSolver::StokesSolver(TopologyHandler *input_topol_handler) : MultiBlockSolver(input_topol_handler), minus_one(-1.0) { @@ -304,7 +326,7 @@ bool StokesSolver::BCExistsOnBdr(const int &global_battr_idx) void StokesSolver::SetupBCOperators() { assert(IsBdrTypeDefined()); - + SetupRHSBCOperators(); SetupDomainBCOperators(); @@ -1451,3 +1473,167 @@ void StokesSolver::ComputeBEIntegral( result += Qvec; } } + +void StokesSolver::SetSubsetComplementaryFlux( + const int N, const int M, const int i0, const int j0, + const Array &subset_bdr_attributes, + const Array &subset_bdrtype +) +{ + // assert N from the N x N array + assert(N * N == numSub); + + // Validate input parameters + assert(i0 >= 0 && j0 >= 0); + assert(i0 + M <= N && j0 + M <= N); + + FiniteElementSpace *ufesm = NULL; + ElementTransformation *eltrans = NULL; + VectorCoefficient *ud = NULL; + VectorCoefficient *usol = NULL; + Mesh *mesh = NULL; + + // initializing complementary flux. + subset_flux::del_u = 0.0; + + // Origin and direction function + Vector *x0 = &(subset_flux::x0); + x0->SetSize(dim); + (*x0) = 0.0; + VectorFunctionCoefficient dir_coeff(dim, subset_flux::dir); + + // Determine the center of domain first. + Vector x1(dim), dx1(dim); + x1 = 0.0; dx1 = 0.0; + double area = 0.0; + ConstantCoefficient one(1.0); + for (int mi = 0; mi < M; mi++) + { + for (int mj = 0; mj < M; mj++) + { + int subset_idx = mi * M + mj; + int orig_idx = (i0 + mi) * N + (j0 + mj); + + mesh = meshes[orig_idx]; + ufesm = ufes[orig_idx]; + + for (int i = 0; i < ufesm -> GetNBE(); i++) + { + const int bdr_attr = mesh->GetBdrAttribute(i); + const int subset_idx = subset_bdr_attributes.Find(bdr_attr); + if (subset_idx < 0) { continue; } + if (subset_bdrtype[subset_idx] != BoundaryType::DIRICHLET) { continue; } + + eltrans = ufesm -> GetBdrElementTransformation (i); + area += ComputeBEIntegral(*ufesm->GetBE(i), *eltrans, one); + ComputeBEIntegral(*ufesm->GetBE(i), *eltrans, dir_coeff, dx1); + x1 += dx1; + } // for (int i = 0; i < ufesm -> GetNBE(); i++) + } // for (int mj = 0; mj < M; mj++) + } // for (int mi = 0; mi < M; mi++) + x1 /= area; + + // set the center of domain for direction function. + (*x0) = x1; + + // Evaluate boundary flux \int u_d \dot n dA. + // If global boundary, evaluate from ud_coeffs. + // If internal subset boundary, evaluate from the solution. + // dirflux is evaluated only on internal boundary. + double bflux = 0.0, dirflux = 0.0; + for (int mi = 0; mi < M; mi++) + { + for (int mj = 0; mj < M; mj++) + { + int subset_idx = mi * M + mj; + int orig_idx = (i0 + mi) * N + (j0 + mj); + + mesh = meshes[orig_idx]; + ufesm = ufes[orig_idx]; + usol = new VectorGridFunctionCoefficient(vels[orig_idx]); + + for (int i = 0; i < ufesm -> GetNBE(); i++) + { + const int bdr_attr = mesh->GetBdrAttribute(i); + const int global_idx = global_bdr_attributes.Find(bdr_attr); + const int subset_idx = subset_bdr_attributes.Find(bdr_attr); + if ((global_idx < 0) && (subset_idx < 0)) continue; + if (subset_bdrtype[subset_idx] != BoundaryType::DIRICHLET) { continue; } + + if ((global_idx >= 0) && (subset_idx >= 0)) // global boundary + ud = ud_coeffs[global_idx]; + else if ((global_idx < 0) && (subset_idx >= 0)) // internal boundary + ud = usol; + else + mfem_error("StokesSolver::SetSubsetComplementaryFlux- subset is not properly created!\n"); + + eltrans = ufesm -> GetBdrElementTransformation (i); + bflux += ComputeBEFlux(*ufesm->GetBE(i), *eltrans, *ud); + + if ((global_idx < 0) && (subset_idx >= 0)) + dirflux += ComputeBEFlux(*ufesm->GetBE(i), *eltrans, dir_coeff); + } // for (int i = 0; i < ufesm -> GetNBE(); i++) + + delete usol; + } // for (int mj = 0; mj < M; mj++) + } // for (int mi = 0; mi < M; mi++) + assert(dirflux > 0.0); + + // Set the flux to ensure incompressibility. + subset_flux::del_u = bflux / dirflux; + VectorFunctionCoefficient subset_flux_coeff(dim, subset_flux::flux); + for (int mi = 0; mi < M; mi++) + { + for (int mj = 0; mj < M; mj++) + { + int subset_idx = mi * M + mj; + int orig_idx = (i0 + mi) * N + (j0 + mj); + + GridFunction tmp(ufes[orig_idx]); + tmp.ProjectCoefficient(subset_flux_coeff); + (*vels[orig_idx]) -= tmp; + } + } + + // Make sure the resulting flux is zero. + double threshold = 1.0e-12; + bflux = 0.0; + for (int mi = 0; mi < M; mi++) + { + for (int mj = 0; mj < M; mj++) + { + int subset_idx = mi * M + mj; + int orig_idx = (i0 + mi) * N + (j0 + mj); + + mesh = meshes[orig_idx]; + ufesm = ufes[orig_idx]; + usol = new VectorGridFunctionCoefficient(vels[orig_idx]); + + for (int i = 0; i < ufesm -> GetNBE(); i++) + { + const int bdr_attr = mesh->GetBdrAttribute(i); + const int global_idx = global_bdr_attributes.Find(bdr_attr); + const int subset_idx = subset_bdr_attributes.Find(bdr_attr); + if ((global_idx < 0) && (subset_idx < 0)) continue; + if (subset_bdrtype[subset_idx] != BoundaryType::DIRICHLET) { continue; } + + if ((global_idx >= 0) && (subset_idx >= 0)) // global boundary + ud = ud_coeffs[global_idx]; + else if ((global_idx < 0) && (subset_idx >= 0)) // internal boundary + ud = usol; + else + mfem_error("StokesSolver::SetSubsetComplementaryFlux- subset is not properly created!\n"); + + eltrans = ufesm -> GetBdrElementTransformation (i); + bflux += ComputeBEFlux(*ufesm->GetBE(i), *eltrans, *ud); + } // for (int i = 0; i < ufesm -> GetNBE(); i++) + + delete usol; + } // for (int mj = 0; mj < M; mj++) + } // for (int mi = 0; mi < M; mi++) + if (abs(bflux) > threshold) + { + printf("boundary flux: %.5E\n", bflux); + mfem_error("StokesSolver::SetSubsetComplementaryFlux- Current boundary setup cannot ensure incompressibility!\n"); + } +} \ No newline at end of file From ff05bb825f1c47a7fda09054e2f23c3fb6e71be7 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 16 Jan 2026 09:27:52 -0800 Subject: [PATCH 16/27] base setup for SchwarzROM. --- CMakeLists.txt | 1 + include/main_workflow.hpp | 3 + include/steady_ns_solver.hpp | 2 + src/main_workflow.cpp | 233 +++++++++++++++++++++++++++++++ src/steady_ns_solver_schwarz.cpp | 134 ++++++++++++++++++ 5 files changed, 373 insertions(+) create mode 100644 src/steady_ns_solver_schwarz.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f9d5c1c5..700b25b3 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/main_workflow.hpp b/include/main_workflow.hpp index c08684b7..381e8af7 100644 --- a/include/main_workflow.hpp +++ b/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); diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 5c35d021..48b841a1 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -198,6 +198,8 @@ friend class SteadyNSOperator; void SaveEQPCoords(const std::string &filename) override; + void SchwarzROM(const int M, const int N, ParameterizedProblem *problem); + private: DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace); diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index a867c061..f054d59b 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -720,6 +720,239 @@ double SingleRun(MPI_Comm comm, const std::string output_file) return error.Max(); } + +double SingleSchwarzRun(MPI_Comm comm, const std::string output_file) +{ + std::string solver_type = config.GetRequiredOption("main/solver"); + if (solver_type != "steady-ns") + mfem_error("SingleSchwarzRun currently only supports steady Navier-Stokes solver!\n"); + + if (config.GetOption("single_run/choose_from_random_sample", false)) + { + RandomSampleGenerator *generator = new RandomSampleGenerator(comm); + generator->SetParamSpaceSizes(); + int idx = UniformRandom(0, generator->GetTotalSampleSize()-1); + // NOTE: this will change config.dict_ + generator->SetSampleParams(idx); + delete generator; + } + + ParameterizedProblem *problem = InitParameterizedProblem(); + SteadyNSSolver *test = new SteadyNSSolver(); + test->InitVariables(); + if (test->UseRom()) test->InitROMHandler(); + test->InitVisualization(); + + StopWatch solveTimer; + std::string solveType = (test->UseRom()) ? "ROM" : "FOM"; + + problem->SetSingleRun(); + assert(test->SetParameterizedProblem(problem)); + + assert(test->IsBdrTypeDefined()); + + int M = config.GetRequiredOption("schwarz/local_size"); + int N = config.GetRequiredOption("schwarz/global_size"); + test->SchwarzROM(M, N, problem); + +// // TODO: there are skippable operations depending on rom/fom mode. +// test->BuildRHSOperators(); +// test->SetupRHSBCOperators(); +// test->AssembleRHS(); + +// const int num_var = test->GetNumVar(); +// Vector rom_assemble(1), rom_solve(1), fom_assemble(1), fom_solve(1), error(num_var); +// rom_assemble = -1.0; rom_solve = -1.0; +// fom_assemble = -1.0; fom_solve = -1.0; +// error = -1.0; + +// ROMHandlerBase *rom = NULL; +// if (test->UseRom()) +// { +// rom = test->GetROMHandler(); +// test->LoadReducedBasis(); + +// if (test->IsNonlinear()) +// test->AllocateROMNlinElems(); +// } + +// solveTimer.Start(); +// if (test->UseRom()) +// { +// printf("ROM with "); +// ROMBuildingLevel save_operator = rom->GetBuildingLevel(); +// TopologyHandlerMode topol_mode = test->GetTopologyMode(); + +// if (topol_mode == TopologyHandlerMode::SUBMESH) +// printf("using SubMesh topology.\n"); +// else if (topol_mode == TopologyHandlerMode::COMPONENT) +// printf("using Component-wise topology.\n"); +// else +// mfem_error("Unknown TopologyHandler Mode!\n"); + +// std::string filename = rom->GetOperatorPrefix() + ".h5"; +// if (save_operator == ROMBuildingLevel::COMPONENT) +// { +// if (topol_mode == TopologyHandlerMode::SUBMESH) +// mfem_error("Submesh does not support component rom building level!\n"); + +// printf("Loading ROM projected elements.. "); +// test->LoadROMLinElems(filename); +// printf("Done!\n"); + +// printf("Assembling ROM linear matrix.. "); +// test->AssembleROMMat(); +// printf("Done!\n"); + +// if (test->IsNonlinear()) +// { +// test->LoadROMNlinElems(rom->GetOperatorPrefix()); +// test->AssembleROMNlinOper(); +// } +// } // if (save_operator == ROMBuildingLevel::COMPONENT) +// else if (save_operator == ROMBuildingLevel::GLOBAL) +// { +// printf("Loading global operator file.. "); +// test->LoadROMOperatorFromFile(filename); +// printf("Done!\n"); +// } // if (save_operator == ROMBuildingLevel::GLOBAL) +// else if (save_operator == ROMBuildingLevel::NONE) +// { +// printf("Building operator file all the way from FOM.. "); +// test->BuildDomainOperators(); +// test->SetupDomainBCOperators(); +// test->AssembleOperator(); +// test->ProjectOperatorOnReducedBasis(); +// printf("Done!\n"); +// } // if (save_operator == ROMBuildingLevel::NONE) +// else +// mfem_error("SingleRun - Unknown ROMBuildingLevel!\n"); + +// printf("Projecting RHS to ROM.. "); +// test->ProjectRHSOnReducedBasis(); +// printf("Done!\n"); +// } // if (test->UseRom()) +// else +// { +// test->BuildDomainOperators(); +// test->SetupDomainBCOperators(); +// test->AssembleOperator(); +// } // not if (test->UseRom()) +// solveTimer.Stop(); +// printf("%s-assemble time: %f seconds.\n", solveType.c_str(), solveTimer.RealTime()); + +// if (test->UseRom()) +// rom_assemble = solveTimer.RealTime(); +// else +// fom_assemble = solveTimer.RealTime(); + +// solveTimer.Clear(); +// solveTimer.Start(); +// if (test->UseRom()) +// { +// test->SolveROM(); +// } +// else +// { +// // TODO: move matrix assembly to here. +// test->Solve(); +// } +// solveTimer.Stop(); +// printf("%s-solve time: %f seconds.\n", solveType.c_str(), solveTimer.RealTime()); + +// if (test->UseRom()) +// rom_solve = solveTimer.RealTime(); +// else +// fom_solve = solveTimer.RealTime(); + +// /* save the ROM system for analysis/debug */ +// bool save_rom = config.GetOption("model_reduction/save_linear_system/enabled", false); +// if (save_rom) +// { +// std::string rom_prefix = config.GetRequiredOption("model_reduction/save_linear_system/prefix"); +// rom->SaveRomSystem(rom_prefix); +// } + +// bool compare_sol = config.GetOption("model_reduction/compare_solution/enabled", false); +// bool load_sol = config.GetOption("model_reduction/compare_solution/load_solution", false); +// if (test->UseRom() && compare_sol) +// { +// BlockVector *romU = test->GetSolutionCopy(); + +// if (load_sol) +// { +// printf("Comparing with the existing FOM solution.\n"); +// std::string fom_file = config.GetRequiredOption("model_reduction/compare_solution/fom_solution_file"); +// test->LoadSolution(fom_file); +// } +// else +// { +// solveTimer.Clear(); +// solveTimer.Start(); +// test->BuildDomainOperators(); +// test->SetupDomainBCOperators(); +// test->AssembleOperator(); +// solveTimer.Stop(); +// printf("FOM-assembly time: %f seconds.\n", solveTimer.RealTime()); +// fom_assemble = solveTimer.RealTime(); + +// solveTimer.Clear(); +// solveTimer.Start(); +// test->Solve(); +// solveTimer.Stop(); +// printf("FOM-solve time: %f seconds.\n", solveTimer.RealTime()); +// fom_solve = solveTimer.RealTime(); +// } + +// test->CompareSolution(*romU, error); + +// bool save_reduced_sol = config.GetOption("model_reduction/compare_solution/save_reduced_solution", false); +// if (save_reduced_sol) +// { +// ROMHandlerBase *rom = test->GetROMHandler(); +// rom->SaveReducedSolution("rom_reduced_sol.txt"); + +// // use ROMHandler::reduced_rhs as a temporary variable. +// rom->ProjectRHSOnReducedBasis(test->GetSolution()); +// rom->SaveReducedRHS("fom_reduced_sol.txt"); +// } + +// // Recover the original ROM solution. +// test->CopySolution(romU); + +// delete romU; +// } + +// // save results to output file. +// if (output_file.length() > 0) +// { +// hid_t file_id; +// herr_t errf = 0; +// file_id = H5Fcreate(output_file.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); +// assert(file_id >= 0); + +// hdf5_utils::WriteDataset(file_id, "rom_assemble", rom_assemble); +// hdf5_utils::WriteDataset(file_id, "rom_solve", rom_solve); +// hdf5_utils::WriteDataset(file_id, "fom_assemble", fom_assemble); +// hdf5_utils::WriteDataset(file_id, "fom_solve", fom_solve); +// hdf5_utils::WriteDataset(file_id, "rel_error", error); + +// errf = H5Fclose(file_id); +// assert(errf >= 0); +// } + +// // Save solution and visualization. +// test->SaveSolution(); +// test->SaveVisualization(); + + delete test; + delete problem; + + // return the maximum error over all variables. + // return error.Max(); + return 0.; +} + void PrintEQPCoords(MPI_Comm comm) { ParameterizedProblem *problem = InitParameterizedProblem(); diff --git a/src/steady_ns_solver_schwarz.cpp b/src/steady_ns_solver_schwarz.cpp new file mode 100644 index 00000000..aad58e23 --- /dev/null +++ b/src/steady_ns_solver_schwarz.cpp @@ -0,0 +1,134 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include "steady_ns_solver.hpp" +#include "component_topology_handler.hpp" +#include "hyperreduction_integ.hpp" +#include "nonlinear_integ.hpp" +// #include "input_parser.hpp" +// #include "hdf5_utils.hpp" +// #include "linalg_utils.hpp" +// #include "dg_bilinear.hpp" +#include "dg_linear.hpp" +#include "etc.hpp" + +using namespace std; +using namespace mfem; + +void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem *problem) +{ + assert(use_rom); + assert(topol_mode == TopologyHandlerMode::COMPONENT); + assert(N * N == numSub); + ComponentTopologyHandler *comp_topol = static_cast(topol_handler); + + int Ns = N - M + 1; + Array sub_topols(Ns * Ns); // will be owned by sub MultiBlockSolvers defined subsequently. + Array sub_solvers(Ns * Ns); + Array *> subset2orig(Ns * Ns); + for (int i = 0; i < Ns; i++) + for (int j = 0; j < Ns; j++) + { + int index = i * Ns + j; + subset2orig[index] = new Array; + sub_topols[index] = new ComponentTopologyHandler(comp_topol, i, j, N, M, *subset2orig[index]); + sub_solvers[index] = new SteadyNSSolver(sub_topols[index]); + } + + // Base initialization and boundary setup by the global problem. + for (int k = 0; k < Ns * Ns; k++) + { + sub_solvers[k]->InitVariables(); + // sub_solvers[k]->InitVisualization(); + sub_solvers[k]->SetParameterizedProblem(problem); + } + + // Setting Dirichlet on internal boundaries. + Array ensure_incomp(Ns * Ns); + ensure_incomp = true; + for (int k = 0; k < Ns * Ns; k++) + { + for (int b = 0; b < sub_solvers[k]->global_bdr_attributes.Size(); b++) + { + if (sub_solvers[k]->bdr_type[b] == BoundaryType::NUM_BDR_TYPE) + sub_solvers[k]->bdr_type[b] = BoundaryType::DIRICHLET; + } + + // save whether incompressibility should be ensured. + for (int b = 0; b < sub_solvers[k]->global_bdr_attributes.Size(); b++) + { + if (sub_solvers[k]->bdr_type[b] == BoundaryType::NEUMANN) + { + ensure_incomp[k] = false; + break; + } + } + } + + // Setting global solution as boundary condition pointer. + Array *> internal_bdr_funcs(Ns * Ns); + Array *> internal_bdr_meshes(Ns * Ns); + Array *> internal_bdr_attr(Ns * Ns); + for (int k = 0; k < Ns * Ns; k++) + { + internal_bdr_funcs[k] = new Array(0); + internal_bdr_meshes[k] = new Array(0); + internal_bdr_attr[k] = new Array(0); + + for (int b = 0; b < sub_solvers[k]->global_bdr_attributes.Size(); b++) + { + // Internal boundary is Dirichlet BC without a function coefficient yet. + bool dirichlet = (sub_solvers[k]->bdr_type[b] == BoundaryType::DIRICHLET); + bool no_func = (!sub_solvers[k]->ud_coeffs[b]); + if (!(dirichlet && no_func)) + continue; + + // Find all the meshes in the subset that have the boundary attribute. + int battr = sub_solvers[k]->global_bdr_attributes[b]; + for (int m = 0; m < M * M; m++) + { + Mesh *mesh = sub_topols[k]->GetMesh(m); + int idx = mesh->bdr_attributes.Find(battr); + if (idx < 0) + continue; + + int global_m = (*subset2orig[k])[m]; + + internal_bdr_meshes[k]->Append(m); + internal_bdr_attr[k]->Append(battr); + internal_bdr_funcs[k]->Append(new VectorGridFunctionCoefficient(vels[global_m])); + } + } + } + + // Assemble ROM operators for sub_solvers. + for (int k = 0; k < Ns * Ns; k++) + { + printf("=== Assembly of %d-th sub_solver ===\n", k+1); + ROMHandlerBase *rom = sub_solvers[k]->GetROMHandler(); + sub_solvers[k]->LoadReducedBasis(); + + // SteadyNS is always nonlinear. + // NOTE: will need if-statement for equation generalization. + sub_solvers[k]->AllocateROMNlinElems(); + + ROMBuildingLevel save_operator = rom->GetBuildingLevel(); + if (save_operator != ROMBuildingLevel::COMPONENT) + mfem_error("SteadyNSSolver::SchwarzROM- SchwarzROM only supports component-level ROM building!\n"); + + printf("Loading ROM projected elements.. "); + sub_solvers[k]->LoadROMLinElems(filename); + printf("Done!\n"); + + printf("Assembling ROM linear matrix.. "); + sub_solvers[k]->AssembleROMMat(); + printf("Done!\n"); + + sub_solvers[k]->LoadROMNlinElems(rom->GetOperatorPrefix()); + sub_solvers[k]->AssembleROMNlinOper(); + } + + DeletePointers(sub_solvers); + DeletePointers(subset2orig); +} \ No newline at end of file From 99baea2f3cee13010b725d686e8db3219e28acdd Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Fri, 16 Jan 2026 10:01:13 -0800 Subject: [PATCH 17/27] minor fix --- bin/main.cpp | 1 + src/component_topology_handler.cpp | 2 +- src/steady_ns_solver_schwarz.cpp | 5 ++++- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/bin/main.cpp b/bin/main.cpp index e84e19a3..e9d57c28 100644 --- a/bin/main.cpp +++ b/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 { diff --git a/src/component_topology_handler.cpp b/src/component_topology_handler.cpp index 2ffd2e89..1d55ace4 100644 --- a/src/component_topology_handler.cpp +++ b/src/component_topology_handler.cpp @@ -1018,7 +1018,7 @@ ComponentTopologyHandler::ComponentTopologyHandler( tf_ptr(global->tf_ptr), inv_tf_ptr(global->inv_tf_ptr) { // assert N from the N x N array - assert(N * N == numSub); + assert(N * N == global->numSub); // Validate input parameters assert(i0 >= 0 && j0 >= 0); diff --git a/src/steady_ns_solver_schwarz.cpp b/src/steady_ns_solver_schwarz.cpp index aad58e23..30c13b71 100644 --- a/src/steady_ns_solver_schwarz.cpp +++ b/src/steady_ns_solver_schwarz.cpp @@ -105,7 +105,9 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * // Assemble ROM operators for sub_solvers. for (int k = 0; k < Ns * Ns; k++) { - printf("=== Assembly of %d-th sub_solver ===\n", k+1); + sub_solvers[k]->InitROMHandler(); + + printf("\n=== Assembly of %d-th sub_solver ===\n\n", k+1); ROMHandlerBase *rom = sub_solvers[k]->GetROMHandler(); sub_solvers[k]->LoadReducedBasis(); @@ -118,6 +120,7 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * mfem_error("SteadyNSSolver::SchwarzROM- SchwarzROM only supports component-level ROM building!\n"); printf("Loading ROM projected elements.. "); + std::string filename = rom->GetOperatorPrefix() + ".h5"; sub_solvers[k]->LoadROMLinElems(filename); printf("Done!\n"); From ed7ca86df7266e7beabb67c858ae2dd29c7513d3 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 16 Jan 2026 12:10:34 -0800 Subject: [PATCH 18/27] schwarz loop. --- include/steady_ns_solver.hpp | 7 +- include/stokes_solver.hpp | 1 + src/main_workflow.cpp | 4 +- src/steady_ns_solver_schwarz.cpp | 111 ++++++++++++++++++++++++++++++- 4 files changed, 120 insertions(+), 3 deletions(-) diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 48b841a1..e4103e89 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -198,7 +198,8 @@ friend class SteadyNSOperator; void SaveEQPCoords(const std::string &filename) override; - void SchwarzROM(const int M, const int N, ParameterizedProblem *problem); + void SchwarzROM(const int M, const int N, ParameterizedProblem *problem, + const int maxIter, const double threshold=1e-15); private: DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace); @@ -212,6 +213,10 @@ friend class SteadyNSOperator; void SaveEQPElems(const std::string &filename); void LoadEQPElems(const std::string &filename); void AssembleROMEQPOper(); + +protected: + void SetupSubsetRHSBCOperators(const Array *meshes, const Array *battrs, + const Array *bfuncs); }; #endif diff --git a/include/stokes_solver.hpp b/include/stokes_solver.hpp index 54743794..f7c966db 100644 --- a/include/stokes_solver.hpp +++ b/include/stokes_solver.hpp @@ -217,6 +217,7 @@ friend class ParameterizedProblem; void ComputeBEIntegral(const FiniteElement &el, ElementTransformation &Tr, VectorCoefficient &Q, Vector &result); +protected: // Assumes meshes are arranged in an N x N square array. // Extracts an M x M subset starting at (i0, j0). // to ensure incompressibility for the subset domain with all velocity dirichlet bc. diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index f054d59b..c4dd734d 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -753,7 +753,9 @@ double SingleSchwarzRun(MPI_Comm comm, const std::string output_file) int M = config.GetRequiredOption("schwarz/local_size"); int N = config.GetRequiredOption("schwarz/global_size"); - test->SchwarzROM(M, N, problem); + int maxIter = config.GetRequiredOption("schwarz/maximum_iteration"); + double threshold = config.GetRequiredOption("schwarz/threshold"); + test->SchwarzROM(M, N, problem, maxIter, threshold); // // TODO: there are skippable operations depending on rom/fom mode. // test->BuildRHSOperators(); diff --git a/src/steady_ns_solver_schwarz.cpp b/src/steady_ns_solver_schwarz.cpp index 30c13b71..1bc71a04 100644 --- a/src/steady_ns_solver_schwarz.cpp +++ b/src/steady_ns_solver_schwarz.cpp @@ -16,7 +16,8 @@ using namespace std; using namespace mfem; -void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem *problem) +void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem *problem, + const int maxIter, const double threshold) { assert(use_rom); assert(topol_mode == TopologyHandlerMode::COMPONENT); @@ -26,11 +27,14 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * int Ns = N - M + 1; Array sub_topols(Ns * Ns); // will be owned by sub MultiBlockSolvers defined subsequently. Array sub_solvers(Ns * Ns); + Array i0s(Ns * Ns), j0s(Ns * Ns); Array *> subset2orig(Ns * Ns); for (int i = 0; i < Ns; i++) for (int j = 0; j < Ns; j++) { int index = i * Ns + j; + i0s[index] = i; + j0s[index] = j; subset2orig[index] = new Array; sub_topols[index] = new ComponentTopologyHandler(comp_topol, i, j, N, M, *subset2orig[index]); sub_solvers[index] = new SteadyNSSolver(sub_topols[index]); @@ -102,6 +106,17 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * } } + // Set up sub_solvers FOM RHS BC operators. + for (int k = 0; k < Ns * Ns; k++) + { + sub_solvers[k]->BuildRHSOperators(); + sub_solvers[k]->SetupRHSBCOperators(); + // Setup RHS BC operator for internal boundaries + sub_solvers[k]->SetupSubsetRHSBCOperators(internal_bdr_meshes[k], + internal_bdr_attr[k], + internal_bdr_funcs[k]); + } + // Assemble ROM operators for sub_solvers. for (int k = 0; k < Ns * Ns; k++) { @@ -131,7 +146,101 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * sub_solvers[k]->LoadROMNlinElems(rom->GetOperatorPrefix()); sub_solvers[k]->AssembleROMNlinOper(); } + + // Global solution initialization (tiny random perturbation) + for (int k = 0; k < U->Size(); k++) + (*U)[k] = 1.0e-5 * UniformRandom(); + + // Main Schwarz loop. + double error = 1.0; + for (int iter = 0; iter < maxIter; iter++) + { + // Sweep through sub-solvers. + for (int k = 0; k < Ns * Ns; k++) + { + // Adjust global solution to ensure divergence-free BC. + if (ensure_incomp[k]) + SetSubsetComplementaryFlux(N, M, i0s[k], j0s[k], + sub_solvers[k]->global_bdr_attributes, + sub_solvers[k]->bdr_type); + + // Assemble FOM-level RHS. + // All RHS BC operators are already defined, + // and linked to the adjusted global solution. + sub_solvers[k]->AssembleRHS(); + + printf("%d-th sub_solver: Projecting RHS to ROM.. ", k+1); + sub_solvers[k]->ProjectRHSOnReducedBasis(); + printf("Done!\n"); + + // Solve for the subsolver + sub_solvers[k]->SolveROM(); + + // Compute relative error after iteration. + double error1 = 0.0; + int norm = 0.0; + for (int m = 0; m < sub_solvers[k]->numSub; m++) + { + double subdomain_error, subdomain_norm; + + const int orig_idx = (*subset2orig[k])[m]; + ComputeSubdomainErrorAndNorm(vels[orig_idx], sub_solvers[k]->vels[m], + subdomain_error, subdomain_norm); + norm += subdomain_norm * subdomain_norm; + error += subdomain_error * subdomain_error; + } + norm = sqrt(norm); + error1 = sqrt(error1); + error1 /= norm; + error = max(error, error1); + + // Project subsolver solution to global solution. + for (int m = 0; m < sub_solvers[k]->numSub; m++) + { + const int orig_idx = (*subset2orig[k])[m]; + (*vels[orig_idx]) = (*(sub_solvers[k]->vels[m])); + } + } // for (int k = 0; k < Ns * Ns; k++) + + // Exit the iterations if error is below threshold. + if (error <= threshold) + { + printf("SteadyNSSolver::SchwarzROM- Schwarz iteration converged. Iteration error: %.4e\n", error); + break; + } + } // for (int iter = 0; iter < maxIter; iter++) + + if (error > threshold) + mfem_error("SteadyNSSolver::SchwarzROM- Schwarz iteration failed to converge!\n"); DeletePointers(sub_solvers); DeletePointers(subset2orig); +} + +void SteadyNSSolver::SetupSubsetRHSBCOperators( + const Array *bmeshes, const Array *battrs, + const Array *bfuncs) +{ + const int N = bmeshes->Size(); + assert((battrs->Size() == N) && (bfuncs->Size() == N)); + + for (int k = 0; k < N; k++) + { + const int m = (*bmeshes)[k]; + const int battr = (*battrs)[k]; + VectorGridFunctionCoefficient *bfunc = (*bfuncs)[k]; + const int bidx = meshes[m]->bdr_attributes.Find(battr); + + assert(fs[m] && gs[m]); + assert(bidx >= 0); + assert(bdr_type[bidx] == BoundaryType::DIRICHLET); + assert(!BCExistsOnBdr(bidx)); // For internal boundary, global ud coefficient is not defined. + + fs[m]->AddBdrFaceIntegrator(new DGVectorDirichletLFIntegrator(*bfunc, *nu_coeff, sigma, kappa), *bdr_markers[bidx]); + + if (full_dg) + gs[m]->AddBdrFaceIntegrator(new DGBoundaryNormalLFIntegrator(*bfunc), *bdr_markers[bidx]); + else + gs[m]->AddBoundaryIntegrator(new DGBoundaryNormalLFIntegrator(*bfunc), *bdr_markers[bidx]); + } } \ No newline at end of file From 29f9462cbf1547c52c5c6703654a5f6a01f6f366 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Sat, 17 Jan 2026 14:23:17 -0800 Subject: [PATCH 19/27] MultiBlockSolver::PrintConfiguration --- include/multiblock_solver.hpp | 2 ++ src/multiblock_solver.cpp | 61 +++++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+) diff --git a/include/multiblock_solver.hpp b/include/multiblock_solver.hpp index c23f88f3..362390b7 100644 --- a/include/multiblock_solver.hpp +++ b/include/multiblock_solver.hpp @@ -271,6 +271,8 @@ friend class ParameterizedProblem; return true; } + void PrintConfiguration() const; + protected: virtual void AssembleROMMat(BlockMatrix &romMat); }; diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index 88ecab56..f02d0ed2 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -783,3 +783,64 @@ void MultiBlockSolver::CompareSolution(BlockVector &test_U, Vector &error) DeletePointers(test_us); } + +void MultiBlockSolver::PrintConfiguration() const +{ + printf("\n\n======= Mesh/Boundary Configuration =======\n\n"); + + printf("--- Global Boundary Attributes ---\n"); + assert(global_bdr_attributes.Size() == numBdr); + assert(bdr_type.Size() == numBdr); + for (int b = 0; b < numBdr; b++) + { + std::string type; + switch (bdr_type[b]) + { + case BoundaryType::ZERO: + { + type = "ZERO"; + break; + } + case BoundaryType::DIRICHLET: + { + type = "DIRICHLET"; + break; + } + case BoundaryType::NEUMANN: + { + type = "NEUMANN"; + break; + } + default: + { + type = "UNDEFINED"; + break; + } + } + printf("Bdr Attr %d: %s\n", global_bdr_attributes[b], type.c_str()); + } + + printf("\n--- Mesh Boundary Attribute Mapping ---\n"); + for (int m = 0; m < numSub; m++) + { + printf("\n------------------------------\n"); + printf(" Mesh %d \n", m); + printf("------------------------------\n"); + Array *bdr_c2g = topol_handler->GetBdrAttrComponentToGlobalMap(m); + printf("Comp:\t"); + for (int k = 0; k < bdr_c2g->Size(); k++) + printf("%d\t", k+1); + printf("\n"); + printf("Glob:\t"); + for (int k = 0; k < bdr_c2g->Size(); k++) + printf("%d\t", (*bdr_c2g)[k]); + printf("\n"); + } + printf("\n"); + + printf("\n--- Port Info ---\n"); + topol_handler->PrintPortInfo(); + + printf("\n\n===========================================\n\n"); + return; +} \ No newline at end of file From a1641469d58453668bd5112d270022fcf0988ede Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Sat, 17 Jan 2026 14:28:09 -0800 Subject: [PATCH 20/27] PrintBoundaryType --- include/parameterized_problem.hpp | 2 ++ src/multiblock_solver.cpp | 25 +------------------------ src/parameterized_problem.cpp | 29 +++++++++++++++++++++++++++++ src/stokes_solver.cpp | 13 +++++-------- 4 files changed, 37 insertions(+), 32 deletions(-) diff --git a/include/parameterized_problem.hpp b/include/parameterized_problem.hpp index 70d287d9..849845a0 100644 --- a/include/parameterized_problem.hpp +++ b/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; diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index f02d0ed2..8d391c51 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -793,30 +793,7 @@ void MultiBlockSolver::PrintConfiguration() const assert(bdr_type.Size() == numBdr); for (int b = 0; b < numBdr; b++) { - std::string type; - switch (bdr_type[b]) - { - case BoundaryType::ZERO: - { - type = "ZERO"; - break; - } - case BoundaryType::DIRICHLET: - { - type = "DIRICHLET"; - break; - } - case BoundaryType::NEUMANN: - { - type = "NEUMANN"; - break; - } - default: - { - type = "UNDEFINED"; - break; - } - } + std::string type = PrintBoundaryType(bdr_type[b]); printf("Bdr Attr %d: %s\n", global_bdr_attributes[b], type.c_str()); } diff --git a/src/parameterized_problem.cpp b/src/parameterized_problem.cpp index 10699fc0..eeb6cbe5 100644 --- a/src/parameterized_problem.cpp +++ b/src/parameterized_problem.cpp @@ -424,6 +424,35 @@ void body_force(const Vector &x, double t, Vector &u) } // namespace function_factory +const std::string PrintBoundaryType(const BoundaryType bdr_type) +{ + std::string type; + switch (bdr_type) + { + case BoundaryType::ZERO: + { + type = "ZERO"; + break; + } + case BoundaryType::DIRICHLET: + { + type = "DIRICHLET"; + break; + } + case BoundaryType::NEUMANN: + { + type = "NEUMANN"; + break; + } + default: + { + type = "UNDEFINED"; + break; + } + } + return type; +} + ParameterizedProblem::ParameterizedProblem() : problem_name(config.GetRequiredOption("parameterized_problem/name")) { diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index 67c19492..f2ea99de 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -1155,15 +1155,12 @@ bool StokesSolver::SetParameterizedProblem(ParameterizedProblem *problem) { Array nz_dbcs(numBdr); nz_dbcs = true; - for (int b = 0; b < problem->battr.Size(); b++) + for (int b = 0; b < global_bdr_attributes.Size(); b++) { - if (problem->bdr_type[b] == BoundaryType::ZERO) - { - if (problem->battr[b] == -1) - { nz_dbcs = false; break; } - else - nz_dbcs[b] = false; - } + if (bdr_type[b] == BoundaryType::ZERO) + nz_dbcs[b] = false; + else + assert(bdr_type[b] == BoundaryType::DIRICHLET); } SetComplementaryFlux(nz_dbcs); } From 38aec96eef57a0b99e89dcf2d90c133fc47d6e26 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Sat, 17 Jan 2026 19:13:46 -0800 Subject: [PATCH 21/27] SteadyNSSolver::SchwarzROM fixed. May need to try a few more things to decrease error level. --- include/steady_ns_solver.hpp | 2 ++ include/stokes_solver.hpp | 3 +- src/parameterized_problem.cpp | 6 ++-- src/steady_ns_solver.cpp | 22 +++++++++++++ src/steady_ns_solver_schwarz.cpp | 56 ++++++++++++++++++++++---------- src/stokes_solver.cpp | 30 +++++++++++++---- 6 files changed, 91 insertions(+), 28 deletions(-) diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index e4103e89..08428824 100644 --- a/include/steady_ns_solver.hpp +++ b/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 diff --git a/include/stokes_solver.hpp b/include/stokes_solver.hpp index f7c966db..670a84c7 100644 --- a/include/stokes_solver.hpp +++ b/include/stokes_solver.hpp @@ -223,7 +223,8 @@ friend class ParameterizedProblem; // to ensure incompressibility for the subset domain with all velocity dirichlet bc. void SetSubsetComplementaryFlux(const int N, const int M, const int i0, const int j0, const Array &subset_bdr_attributes, - const Array &subset_bdrtype); + const Array &subset_bdrtype, + const ParameterizedProblem *problem); }; #endif diff --git a/src/parameterized_problem.cpp b/src/parameterized_problem.cpp index eeb6cbe5..f8896465 100644 --- a/src/parameterized_problem.cpp +++ b/src/parameterized_problem.cpp @@ -148,9 +148,9 @@ void ubdr(const Vector &x, double t, Vector &y) } // ensure incompressibility. - Vector del_u(dim); - flow_problem::flux(x, del_u); - y -= del_u; + Vector del_flux(dim); + flow_problem::flux(x, del_flux); + y -= del_flux; } } // namespace component_flow diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index 560bf8bd..413afa16 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -241,6 +241,28 @@ Operator& SteadyNSTensorROM::GetGradient(const Vector &x) const return *jac_mono; } +void SteadyNSTensorROM::SaveOperator(std::string filename) +{ + hid_t file_id; + herr_t errf = 0; + file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + assert(file_id >= 0); + + DenseMatrix linOp; + linearOp->ToDenseMatrix(linOp); + hdf5_utils::WriteDataset(file_id, "linear_op", linOp); + + for (int m = 0; m < numSub; m++) + { + std::string dsetname = "tensor" + std::to_string(m); + hdf5_utils::WriteDataset(file_id, dsetname.c_str(), *hs[m]); + } + + errf = H5Fclose(file_id); + assert(errf >= 0); + return; +} + /* SteadyNSTensorROM */ diff --git a/src/steady_ns_solver_schwarz.cpp b/src/steady_ns_solver_schwarz.cpp index 1bc71a04..637642e6 100644 --- a/src/steady_ns_solver_schwarz.cpp +++ b/src/steady_ns_solver_schwarz.cpp @@ -6,10 +6,6 @@ #include "component_topology_handler.hpp" #include "hyperreduction_integ.hpp" #include "nonlinear_integ.hpp" -// #include "input_parser.hpp" -// #include "hdf5_utils.hpp" -// #include "linalg_utils.hpp" -// #include "dg_bilinear.hpp" #include "dg_linear.hpp" #include "etc.hpp" @@ -84,6 +80,7 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * { // Internal boundary is Dirichlet BC without a function coefficient yet. bool dirichlet = (sub_solvers[k]->bdr_type[b] == BoundaryType::DIRICHLET); + dirichlet = dirichlet || (sub_solvers[k]->bdr_type[b] == BoundaryType::ZERO); bool no_func = (!sub_solvers[k]->ud_coeffs[b]); if (!(dirichlet && no_func)) continue; @@ -147,31 +144,54 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * sub_solvers[k]->AssembleROMNlinOper(); } - // Global solution initialization (tiny random perturbation) - for (int k = 0; k < U->Size(); k++) - (*U)[k] = 1.0e-5 * UniformRandom(); + // Global solution initialization + // HACK: we assume the ud_coeff is the same for all non-zero Dirichlet condition. + for (int b = 0; b < global_bdr_attributes.Size(); b++) + { + if ((bdr_type[b] == BoundaryType::DIRICHLET) && ud_coeffs[b]) + { + for (int m = 0; m < numSub; m++) + vels[m]->ProjectCoefficient(*ud_coeffs[b]); + break; + } + } + // for (int k = 0; k < U->Size(); k++) + // (*U)[k] = 1.0e-5 * UniformRandom(); // Main Schwarz loop. - double error = 1.0; + bool use_restart = config.GetOption("rom_solver/use_restart", false); + double error = 0.0; for (int iter = 0; iter < maxIter; iter++) { + error = 0.0; // Sweep through sub-solvers. for (int k = 0; k < Ns * Ns; k++) { + // int k = sweep_index[s]; + printf("Switched to %dth subsolver.\n", k+1); + // Adjust global solution to ensure divergence-free BC. if (ensure_incomp[k]) SetSubsetComplementaryFlux(N, M, i0s[k], j0s[k], sub_solvers[k]->global_bdr_attributes, - sub_solvers[k]->bdr_type); + sub_solvers[k]->bdr_type, problem); + + // Project global solution to subsolver solution. + if (use_restart) + { + for (int m = 0; m < sub_solvers[k]->numSub; m++) + { + const int orig_idx = (*subset2orig[k])[m]; + (*(sub_solvers[k]->vels[m])) = (*vels[orig_idx]); + } + } // Assemble FOM-level RHS. // All RHS BC operators are already defined, // and linked to the adjusted global solution. sub_solvers[k]->AssembleRHS(); - printf("%d-th sub_solver: Projecting RHS to ROM.. ", k+1); sub_solvers[k]->ProjectRHSOnReducedBasis(); - printf("Done!\n"); // Solve for the subsolver sub_solvers[k]->SolveROM(); @@ -202,10 +222,11 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * } } // for (int k = 0; k < Ns * Ns; k++) + printf("Iteration %d error: %.4e\n", iter+1, error); // Exit the iterations if error is below threshold. if (error <= threshold) { - printf("SteadyNSSolver::SchwarzROM- Schwarz iteration converged. Iteration error: %.4e\n", error); + printf("SteadyNSSolver::SchwarzROM- Schwarz iteration converged.\n"); break; } } // for (int iter = 0; iter < maxIter; iter++) @@ -229,18 +250,19 @@ void SteadyNSSolver::SetupSubsetRHSBCOperators( const int m = (*bmeshes)[k]; const int battr = (*battrs)[k]; VectorGridFunctionCoefficient *bfunc = (*bfuncs)[k]; + const int global_idx = global_bdr_attributes.Find(battr); const int bidx = meshes[m]->bdr_attributes.Find(battr); assert(fs[m] && gs[m]); assert(bidx >= 0); - assert(bdr_type[bidx] == BoundaryType::DIRICHLET); - assert(!BCExistsOnBdr(bidx)); // For internal boundary, global ud coefficient is not defined. + assert(bdr_type[global_idx] == BoundaryType::DIRICHLET); + assert(!BCExistsOnBdr(global_idx)); // For internal boundary, global ud coefficient is not defined. - fs[m]->AddBdrFaceIntegrator(new DGVectorDirichletLFIntegrator(*bfunc, *nu_coeff, sigma, kappa), *bdr_markers[bidx]); + fs[m]->AddBdrFaceIntegrator(new DGVectorDirichletLFIntegrator(*bfunc, *nu_coeff, sigma, kappa), *bdr_markers[global_idx]); if (full_dg) - gs[m]->AddBdrFaceIntegrator(new DGBoundaryNormalLFIntegrator(*bfunc), *bdr_markers[bidx]); + gs[m]->AddBdrFaceIntegrator(new DGBoundaryNormalLFIntegrator(*bfunc), *bdr_markers[global_idx]); else - gs[m]->AddBoundaryIntegrator(new DGBoundaryNormalLFIntegrator(*bfunc), *bdr_markers[bidx]); + gs[m]->AddBoundaryIntegrator(new DGBoundaryNormalLFIntegrator(*bfunc), *bdr_markers[global_idx]); } } \ No newline at end of file diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index f2ea99de..dabc4021 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -17,6 +17,7 @@ namespace subset_flux { double del_u; Vector x0; + // Vector dir; void dir(const Vector &x, Vector &y) { @@ -30,6 +31,7 @@ namespace subset_flux void flux(const Vector &x, Vector &y) { + // y = dir; dir(x, y); y *= del_u; } @@ -1474,7 +1476,8 @@ void StokesSolver::ComputeBEIntegral( void StokesSolver::SetSubsetComplementaryFlux( const int N, const int M, const int i0, const int j0, const Array &subset_bdr_attributes, - const Array &subset_bdrtype + const Array &subset_bdrtype, + const ParameterizedProblem *problem ) { // assert N from the N x N array @@ -1493,11 +1496,24 @@ void StokesSolver::SetSubsetComplementaryFlux( // initializing complementary flux. subset_flux::del_u = 0.0; - // Origin and direction function - Vector *x0 = &(subset_flux::x0); - x0->SetSize(dim); - (*x0) = 0.0; + // // Determine direction + // // HACK: We assume the first vector bdr ptr of the problem indicates non-zero Dirichlet condition. + // // HACK: We also assume the vector bdr ptr is constant in space. + // { + // subset_flux::dir.SetSize(dim); + + // Vector xtmp(dim); + // (problem->vector_bdr_ptr[0])(xtmp, 0.0, subset_flux::dir); + + // // Normalize + // subset_flux::dir /= subset_flux::dir.Norml2(); + // } + + // direction function + subset_flux::x0.SetSize(dim); + subset_flux::x0 = 0.0; VectorFunctionCoefficient dir_coeff(dim, subset_flux::dir); + VectorFunctionCoefficient subset_flux_coeff(dim, subset_flux::flux); // Determine the center of domain first. Vector x1(dim), dx1(dim); @@ -1531,7 +1547,7 @@ void StokesSolver::SetSubsetComplementaryFlux( x1 /= area; // set the center of domain for direction function. - (*x0) = x1; + subset_flux::x0 = x1; // Evaluate boundary flux \int u_d \dot n dA. // If global boundary, evaluate from ud_coeffs. @@ -1574,11 +1590,11 @@ void StokesSolver::SetSubsetComplementaryFlux( delete usol; } // for (int mj = 0; mj < M; mj++) } // for (int mi = 0; mi < M; mi++) +printf("bflux: %.5e\n", bflux); assert(dirflux > 0.0); // Set the flux to ensure incompressibility. subset_flux::del_u = bflux / dirflux; - VectorFunctionCoefficient subset_flux_coeff(dim, subset_flux::flux); for (int mi = 0; mi < M; mi++) { for (int mj = 0; mj < M; mj++) From 3fe08472b3f91a0284bc8eae737fac2a2aaf3c54 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Sun, 18 Jan 2026 19:45:01 -0800 Subject: [PATCH 22/27] rom2fom error comparison and finalizing the loop.; --- include/steady_ns_solver.hpp | 4 +- src/main_workflow.cpp | 214 ++++++++----------------------- src/steady_ns_solver_schwarz.cpp | 21 ++- 3 files changed, 74 insertions(+), 165 deletions(-) diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 08428824..1433f0a1 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -201,7 +201,9 @@ friend class SteadyNSOperator; void SaveEQPCoords(const std::string &filename) override; void SchwarzROM(const int M, const int N, ParameterizedProblem *problem, - const int maxIter, const double threshold=1e-15); + const int maxIter, const double threshold=1e-15, + const int hist_track=3, const double plateau_range=1e-1, + const bool use_restart=false); private: DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace); diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index c4dd734d..0519f780 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -745,6 +745,8 @@ double SingleSchwarzRun(MPI_Comm comm, const std::string output_file) StopWatch solveTimer; std::string solveType = (test->UseRom()) ? "ROM" : "FOM"; + const int num_var = test->GetNumVar(); + Vector error(num_var); problem->SetSingleRun(); assert(test->SetParameterizedProblem(problem)); @@ -755,175 +757,65 @@ double SingleSchwarzRun(MPI_Comm comm, const std::string output_file) int N = config.GetRequiredOption("schwarz/global_size"); int maxIter = config.GetRequiredOption("schwarz/maximum_iteration"); double threshold = config.GetRequiredOption("schwarz/threshold"); - test->SchwarzROM(M, N, problem, maxIter, threshold); + int hist_track = config.GetOption("schwarz/history_track", 3); + double plateau_range = config.GetOption("schwarz/plateau_range", 1e-1); + bool use_restart = config.GetOption("rom_solver/use_restart", false); + test->SchwarzROM(M, N, problem, maxIter, threshold, + hist_track, plateau_range, use_restart); -// // TODO: there are skippable operations depending on rom/fom mode. -// test->BuildRHSOperators(); -// test->SetupRHSBCOperators(); -// test->AssembleRHS(); + bool compare_sol = config.GetOption("model_reduction/compare_solution/enabled", false); + bool load_sol = config.GetOption("model_reduction/compare_solution/load_solution", false); + if (compare_sol) + { + BlockVector *romU = test->GetSolutionCopy(); -// const int num_var = test->GetNumVar(); -// Vector rom_assemble(1), rom_solve(1), fom_assemble(1), fom_solve(1), error(num_var); -// rom_assemble = -1.0; rom_solve = -1.0; -// fom_assemble = -1.0; fom_solve = -1.0; -// error = -1.0; + if (load_sol) + { + printf("Comparing with the existing FOM solution.\n"); + std::string fom_file = config.GetRequiredOption("model_reduction/compare_solution/fom_solution_file"); + test->LoadSolution(fom_file); + } + else + { + solveTimer.Clear(); + solveTimer.Start(); + test->BuildRHSOperators(); + test->SetupRHSBCOperators(); + test->AssembleRHS(); -// ROMHandlerBase *rom = NULL; -// if (test->UseRom()) -// { -// rom = test->GetROMHandler(); -// test->LoadReducedBasis(); + test->BuildDomainOperators(); + test->SetupDomainBCOperators(); + test->AssembleOperator(); + solveTimer.Stop(); + printf("FOM-assembly time: %f seconds.\n", solveTimer.RealTime()); + // fom_assemble = solveTimer.RealTime(); -// if (test->IsNonlinear()) -// test->AllocateROMNlinElems(); -// } + solveTimer.Clear(); + solveTimer.Start(); + test->Solve(); + solveTimer.Stop(); + printf("FOM-solve time: %f seconds.\n", solveTimer.RealTime()); + // fom_solve = solveTimer.RealTime(); + } -// solveTimer.Start(); -// if (test->UseRom()) -// { -// printf("ROM with "); -// ROMBuildingLevel save_operator = rom->GetBuildingLevel(); -// TopologyHandlerMode topol_mode = test->GetTopologyMode(); - -// if (topol_mode == TopologyHandlerMode::SUBMESH) -// printf("using SubMesh topology.\n"); -// else if (topol_mode == TopologyHandlerMode::COMPONENT) -// printf("using Component-wise topology.\n"); -// else -// mfem_error("Unknown TopologyHandler Mode!\n"); - -// std::string filename = rom->GetOperatorPrefix() + ".h5"; -// if (save_operator == ROMBuildingLevel::COMPONENT) -// { -// if (topol_mode == TopologyHandlerMode::SUBMESH) -// mfem_error("Submesh does not support component rom building level!\n"); - -// printf("Loading ROM projected elements.. "); -// test->LoadROMLinElems(filename); -// printf("Done!\n"); - -// printf("Assembling ROM linear matrix.. "); -// test->AssembleROMMat(); -// printf("Done!\n"); - -// if (test->IsNonlinear()) -// { -// test->LoadROMNlinElems(rom->GetOperatorPrefix()); -// test->AssembleROMNlinOper(); -// } -// } // if (save_operator == ROMBuildingLevel::COMPONENT) -// else if (save_operator == ROMBuildingLevel::GLOBAL) -// { -// printf("Loading global operator file.. "); -// test->LoadROMOperatorFromFile(filename); -// printf("Done!\n"); -// } // if (save_operator == ROMBuildingLevel::GLOBAL) -// else if (save_operator == ROMBuildingLevel::NONE) -// { -// printf("Building operator file all the way from FOM.. "); -// test->BuildDomainOperators(); -// test->SetupDomainBCOperators(); -// test->AssembleOperator(); -// test->ProjectOperatorOnReducedBasis(); -// printf("Done!\n"); -// } // if (save_operator == ROMBuildingLevel::NONE) -// else -// mfem_error("SingleRun - Unknown ROMBuildingLevel!\n"); - -// printf("Projecting RHS to ROM.. "); -// test->ProjectRHSOnReducedBasis(); -// printf("Done!\n"); -// } // if (test->UseRom()) -// else -// { -// test->BuildDomainOperators(); -// test->SetupDomainBCOperators(); -// test->AssembleOperator(); -// } // not if (test->UseRom()) -// solveTimer.Stop(); -// printf("%s-assemble time: %f seconds.\n", solveType.c_str(), solveTimer.RealTime()); - -// if (test->UseRom()) -// rom_assemble = solveTimer.RealTime(); -// else -// fom_assemble = solveTimer.RealTime(); - -// solveTimer.Clear(); -// solveTimer.Start(); -// if (test->UseRom()) -// { -// test->SolveROM(); -// } -// else -// { -// // TODO: move matrix assembly to here. -// test->Solve(); -// } -// solveTimer.Stop(); -// printf("%s-solve time: %f seconds.\n", solveType.c_str(), solveTimer.RealTime()); + test->CompareSolution(*romU, error); + + bool save_reduced_sol = config.GetOption("model_reduction/compare_solution/save_reduced_solution", false); + if (save_reduced_sol) + { + ROMHandlerBase *rom = test->GetROMHandler(); + rom->SaveReducedSolution("rom_reduced_sol.txt"); -// if (test->UseRom()) -// rom_solve = solveTimer.RealTime(); -// else -// fom_solve = solveTimer.RealTime(); + // use ROMHandler::reduced_rhs as a temporary variable. + rom->ProjectRHSOnReducedBasis(test->GetSolution()); + rom->SaveReducedRHS("fom_reduced_sol.txt"); + } -// /* save the ROM system for analysis/debug */ -// bool save_rom = config.GetOption("model_reduction/save_linear_system/enabled", false); -// if (save_rom) -// { -// std::string rom_prefix = config.GetRequiredOption("model_reduction/save_linear_system/prefix"); -// rom->SaveRomSystem(rom_prefix); -// } + // Recover the original ROM solution. + test->CopySolution(romU); -// bool compare_sol = config.GetOption("model_reduction/compare_solution/enabled", false); -// bool load_sol = config.GetOption("model_reduction/compare_solution/load_solution", false); -// if (test->UseRom() && compare_sol) -// { -// BlockVector *romU = test->GetSolutionCopy(); - -// if (load_sol) -// { -// printf("Comparing with the existing FOM solution.\n"); -// std::string fom_file = config.GetRequiredOption("model_reduction/compare_solution/fom_solution_file"); -// test->LoadSolution(fom_file); -// } -// else -// { -// solveTimer.Clear(); -// solveTimer.Start(); -// test->BuildDomainOperators(); -// test->SetupDomainBCOperators(); -// test->AssembleOperator(); -// solveTimer.Stop(); -// printf("FOM-assembly time: %f seconds.\n", solveTimer.RealTime()); -// fom_assemble = solveTimer.RealTime(); - -// solveTimer.Clear(); -// solveTimer.Start(); -// test->Solve(); -// solveTimer.Stop(); -// printf("FOM-solve time: %f seconds.\n", solveTimer.RealTime()); -// fom_solve = solveTimer.RealTime(); -// } - -// test->CompareSolution(*romU, error); - -// bool save_reduced_sol = config.GetOption("model_reduction/compare_solution/save_reduced_solution", false); -// if (save_reduced_sol) -// { -// ROMHandlerBase *rom = test->GetROMHandler(); -// rom->SaveReducedSolution("rom_reduced_sol.txt"); - -// // use ROMHandler::reduced_rhs as a temporary variable. -// rom->ProjectRHSOnReducedBasis(test->GetSolution()); -// rom->SaveReducedRHS("fom_reduced_sol.txt"); -// } - -// // Recover the original ROM solution. -// test->CopySolution(romU); - -// delete romU; -// } + delete romU; + } // // save results to output file. // if (output_file.length() > 0) diff --git a/src/steady_ns_solver_schwarz.cpp b/src/steady_ns_solver_schwarz.cpp index 637642e6..ec06e4fa 100644 --- a/src/steady_ns_solver_schwarz.cpp +++ b/src/steady_ns_solver_schwarz.cpp @@ -13,7 +13,9 @@ using namespace std; using namespace mfem; void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem *problem, - const int maxIter, const double threshold) + const int maxIter, const double threshold, + const int hist_track, const double plateau_range, + const bool use_restart) { assert(use_rom); assert(topol_mode == TopologyHandlerMode::COMPONENT); @@ -159,8 +161,9 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * // (*U)[k] = 1.0e-5 * UniformRandom(); // Main Schwarz loop. - bool use_restart = config.GetOption("rom_solver/use_restart", false); double error = 0.0; + double max_error, min_error; + Array error_hist(0); for (int iter = 0; iter < maxIter; iter++) { error = 0.0; @@ -229,9 +232,21 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * printf("SteadyNSSolver::SchwarzROM- Schwarz iteration converged.\n"); break; } + + error_hist.Prepend(error); + if (error_hist.Size() > hist_track) + { + error_hist.DeleteLast(); + max_error = error_hist.Max(); + min_error = error_hist.Min(); + if ((max_error - min_error) < plateau_range * max_error) + break; + } } // for (int iter = 0; iter < maxIter; iter++) - if (error > threshold) + if ((max_error - min_error) < plateau_range * max_error) + printf("SteadyNSSolver::SchwarzROM- error is plateaued. Exiting the iterations.\n"); + else if (error > threshold) mfem_error("SteadyNSSolver::SchwarzROM- Schwarz iteration failed to converge!\n"); DeletePointers(sub_solvers); From 55dd99d34c40bd9e840c386f15c744cd17caa090 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Sun, 18 Jan 2026 20:30:38 -0800 Subject: [PATCH 23/27] output setup for SchwarzROM. --- include/steady_ns_solver.hpp | 3 +- src/main_workflow.cpp | 65 ++++++++++++++++++++------------ src/steady_ns_solver_schwarz.cpp | 29 +++++++++----- 3 files changed, 62 insertions(+), 35 deletions(-) diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 1433f0a1..16a05f34 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -201,8 +201,9 @@ 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 &error_hist, const int maxIter, const double threshold=1e-15, - const int hist_track=3, const double plateau_range=1e-1, + const int plateau_track=3, const double plateau_range=1e-1, const bool use_restart=false); private: diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index 0519f780..a360a23e 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -753,18 +753,26 @@ double SingleSchwarzRun(MPI_Comm comm, const std::string output_file) assert(test->IsBdrTypeDefined()); + // Schwarz ROM inputs int M = config.GetRequiredOption("schwarz/local_size"); int N = config.GetRequiredOption("schwarz/global_size"); int maxIter = config.GetRequiredOption("schwarz/maximum_iteration"); double threshold = config.GetRequiredOption("schwarz/threshold"); - int hist_track = config.GetOption("schwarz/history_track", 3); + int plateau_track = config.GetOption("schwarz/plateau_track", 3); double plateau_range = config.GetOption("schwarz/plateau_range", 1e-1); bool use_restart = config.GetOption("rom_solver/use_restart", false); - test->SchwarzROM(M, N, problem, maxIter, threshold, - hist_track, plateau_range, use_restart); + // Schwarz ROM outputs + int num_solve = -1; + double rom_solve = -1.0; + Array error_hist(0); + test->SchwarzROM(M, N, problem, rom_solve, num_solve, error_hist, + maxIter, threshold, plateau_track, plateau_range, use_restart); + printf("SchwarzROM solve time: %f seconds.\n", rom_solve); + printf("SchwarzROM number of solve: %d.\n", num_solve); bool compare_sol = config.GetOption("model_reduction/compare_solution/enabled", false); bool load_sol = config.GetOption("model_reduction/compare_solution/load_solution", false); + double fom_solve = -1.0; if (compare_sol) { BlockVector *romU = test->GetSolutionCopy(); @@ -795,7 +803,7 @@ double SingleSchwarzRun(MPI_Comm comm, const std::string output_file) test->Solve(); solveTimer.Stop(); printf("FOM-solve time: %f seconds.\n", solveTimer.RealTime()); - // fom_solve = solveTimer.RealTime(); + fom_solve = solveTimer.RealTime(); } test->CompareSolution(*romU, error); @@ -817,27 +825,34 @@ double SingleSchwarzRun(MPI_Comm comm, const std::string output_file) delete romU; } -// // save results to output file. -// if (output_file.length() > 0) -// { -// hid_t file_id; -// herr_t errf = 0; -// file_id = H5Fcreate(output_file.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); -// assert(file_id >= 0); - -// hdf5_utils::WriteDataset(file_id, "rom_assemble", rom_assemble); -// hdf5_utils::WriteDataset(file_id, "rom_solve", rom_solve); -// hdf5_utils::WriteDataset(file_id, "fom_assemble", fom_assemble); -// hdf5_utils::WriteDataset(file_id, "fom_solve", fom_solve); -// hdf5_utils::WriteDataset(file_id, "rel_error", error); - -// errf = H5Fclose(file_id); -// assert(errf >= 0); -// } - -// // Save solution and visualization. -// test->SaveSolution(); -// test->SaveVisualization(); + // save results to output file. + if (output_file.length() > 0) + { + hid_t file_id; + herr_t errf = 0; + file_id = H5Fcreate(output_file.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + assert(file_id >= 0); + + Array num_solve_output(1); + num_solve_output = num_solve; + Vector rom_solve_output(1), fom_solve_output(1); + rom_solve_output = rom_solve; + fom_solve_output = fom_solve; + + hdf5_utils::WriteDataset(file_id, "rom_num_solve", num_solve_output); + hdf5_utils::WriteDataset(file_id, "rom_solve", rom_solve_output); + // hdf5_utils::WriteDataset(file_id, "fom_assemble", fom_assemble); + hdf5_utils::WriteDataset(file_id, "fom_solve", fom_solve_output); + hdf5_utils::WriteDataset(file_id, "rel_error", error); + hdf5_utils::WriteDataset(file_id, "error_hist", error_hist); + + errf = H5Fclose(file_id); + assert(errf >= 0); + } + + // Save solution and visualization. + test->SaveSolution(); + test->SaveVisualization(); delete test; delete problem; diff --git a/src/steady_ns_solver_schwarz.cpp b/src/steady_ns_solver_schwarz.cpp index ec06e4fa..b30dd9ba 100644 --- a/src/steady_ns_solver_schwarz.cpp +++ b/src/steady_ns_solver_schwarz.cpp @@ -13,13 +13,17 @@ using namespace std; using namespace mfem; void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem *problem, + double &solve_time, int &num_solve, Array &error_hist, const int maxIter, const double threshold, - const int hist_track, const double plateau_range, + const int plateau_track, const double plateau_range, const bool use_restart) { assert(use_rom); assert(topol_mode == TopologyHandlerMode::COMPONENT); assert(N * N == numSub); + StopWatch solveTimer; + solveTimer.Clear(); + error_hist.DeleteAll(); ComponentTopologyHandler *comp_topol = static_cast(topol_handler); int Ns = N - M + 1; @@ -161,9 +165,10 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * // (*U)[k] = 1.0e-5 * UniformRandom(); // Main Schwarz loop. + num_solve = 0; double error = 0.0; double max_error, min_error; - Array error_hist(0); + Array plateau_hist(0); for (int iter = 0; iter < maxIter; iter++) { error = 0.0; @@ -197,7 +202,10 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * sub_solvers[k]->ProjectRHSOnReducedBasis(); // Solve for the subsolver + solveTimer.Start(); sub_solvers[k]->SolveROM(); + solveTimer.Stop(); + num_solve++; // Compute relative error after iteration. double error1 = 0.0; @@ -207,15 +215,16 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * double subdomain_error, subdomain_norm; const int orig_idx = (*subset2orig[k])[m]; - ComputeSubdomainErrorAndNorm(vels[orig_idx], sub_solvers[k]->vels[m], + ComputeSubdomainErrorAndNorm(sub_solvers[k]->vels[m], vels[orig_idx], subdomain_error, subdomain_norm); norm += subdomain_norm * subdomain_norm; - error += subdomain_error * subdomain_error; + error1 += subdomain_error * subdomain_error; } norm = sqrt(norm); error1 = sqrt(error1); error1 /= norm; error = max(error, error1); + error_hist.Append(error1); // Project subsolver solution to global solution. for (int m = 0; m < sub_solvers[k]->numSub; m++) @@ -233,17 +242,19 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * break; } - error_hist.Prepend(error); - if (error_hist.Size() > hist_track) + plateau_hist.Prepend(error); + if (plateau_hist.Size() > plateau_track) { - error_hist.DeleteLast(); - max_error = error_hist.Max(); - min_error = error_hist.Min(); + plateau_hist.DeleteLast(); + max_error = plateau_hist.Max(); + min_error = plateau_hist.Min(); if ((max_error - min_error) < plateau_range * max_error) break; } } // for (int iter = 0; iter < maxIter; iter++) + solve_time = solveTimer.RealTime(); + if ((max_error - min_error) < plateau_range * max_error) printf("SteadyNSSolver::SchwarzROM- error is plateaued. Exiting the iterations.\n"); else if (error > threshold) From 2e45f4672324e3bbc57f4df7c933255f9dd6ee3e Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Sun, 18 Jan 2026 20:43:15 -0800 Subject: [PATCH 24/27] minor fix on MultiBlockSolver::SetParameterizedProblem --- src/multiblock_solver.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index 8d391c51..adec0e12 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -502,6 +502,13 @@ void MultiBlockSolver::SaveVisualization(const int step, const double time) bool MultiBlockSolver::SetParameterizedProblem(ParameterizedProblem *problem) { assert(bdr_type.Size() == global_bdr_attributes.Size()); + // Check if boundary type is set for all boundaries (battr == -1) + if ((problem->battr.Size() == 1) && (problem->battr[0] == -1)) + { + bdr_type = problem->bdr_type[0]; + return true; + } + for (int b = 0; b < global_bdr_attributes.Size(); b++) { int idx = problem->battr.Find(global_bdr_attributes[b]); From 8a16071d9412a704d1702892ff9ac52e1e5d2229 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Mon, 19 Jan 2026 12:01:41 -0800 Subject: [PATCH 25/27] SteadyNSSolver::SchwarzROM option for initial tolerance relaxation. --- include/steady_ns_solver.hpp | 2 +- src/main_workflow.cpp | 3 ++- src/rom_handler.cpp | 3 +++ src/steady_ns_solver_schwarz.cpp | 18 +++++++++++++++++- 4 files changed, 23 insertions(+), 3 deletions(-) diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 16a05f34..2faf8c9f 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -204,7 +204,7 @@ friend class SteadyNSOperator; double &solve_time, int &num_solve, Array &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 bool use_restart=false, const double initial_tol=-1.0); private: DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace); diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index a360a23e..fbc99d7a 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -761,12 +761,13 @@ double SingleSchwarzRun(MPI_Comm comm, const std::string output_file) int plateau_track = config.GetOption("schwarz/plateau_track", 3); double plateau_range = config.GetOption("schwarz/plateau_range", 1e-1); bool use_restart = config.GetOption("rom_solver/use_restart", false); + double initial_tol = config.GetOption("schwarz/initial_tolerance", -1.0); // Schwarz ROM outputs int num_solve = -1; double rom_solve = -1.0; Array error_hist(0); test->SchwarzROM(M, N, problem, rom_solve, num_solve, error_hist, - maxIter, threshold, plateau_track, plateau_range, use_restart); + maxIter, threshold, plateau_track, plateau_range, use_restart, initial_tol); printf("SchwarzROM solve time: %f seconds.\n", rom_solve); printf("SchwarzROM number of solve: %d.\n", num_solve); diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index e86f580d..0e5d0417 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -691,6 +691,9 @@ void MFEMROMHandler::NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec newton_solver.SetMaxIter(maxIter); newton_solver.Mult(*reduced_rhs, *reduced_sol); + bool converged = newton_solver.GetConverged(); + if (!converged) + mfem_error("MFEMROMHandler::NonlinearSolve- Newton solver not converged!\n"); } else if (nlin_solver == "cg") { diff --git a/src/steady_ns_solver_schwarz.cpp b/src/steady_ns_solver_schwarz.cpp index b30dd9ba..ed946774 100644 --- a/src/steady_ns_solver_schwarz.cpp +++ b/src/steady_ns_solver_schwarz.cpp @@ -16,7 +16,7 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * double &solve_time, int &num_solve, Array &error_hist, const int maxIter, const double threshold, const int plateau_track, const double plateau_range, - const bool use_restart) + const bool use_restart, const double initial_tol) { assert(use_rom); assert(topol_mode == TopologyHandlerMode::COMPONENT); @@ -169,8 +169,17 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * double error = 0.0; double max_error, min_error; Array plateau_hist(0); + double tolerance = -1.0; for (int iter = 0; iter < maxIter; iter++) { + // For initial iteration, relax the tolerance for convergence. + if ((iter == 0) && (initial_tol > 0.0)) + { + tolerance = config.GetOption("rom_solver/relative_tolerance", 1.e-10); + config.SetOption("rom_solver/relative_tolerance", initial_tol); + config.SetOption("rom_solver/absolute_tolerance", initial_tol); + } + error = 0.0; // Sweep through sub-solvers. for (int k = 0; k < Ns * Ns; k++) @@ -251,6 +260,13 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * if ((max_error - min_error) < plateau_range * max_error) break; } + + // Recover the tolerance after the initial iteration. + if ((iter == 0) && (initial_tol > 0.0)) + { + config.SetOption("rom_solver/relative_tolerance", tolerance); + config.SetOption("rom_solver/absolute_tolerance", tolerance); + } } // for (int iter = 0; iter < maxIter; iter++) solve_time = solveTimer.RealTime(); From 04c2bfb6a860e9683c97109de86e7126188ea555 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Mon, 19 Jan 2026 13:35:15 -0800 Subject: [PATCH 26/27] minor bug fix. --- src/steady_ns_solver_schwarz.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/steady_ns_solver_schwarz.cpp b/src/steady_ns_solver_schwarz.cpp index ed946774..2af0bdf2 100644 --- a/src/steady_ns_solver_schwarz.cpp +++ b/src/steady_ns_solver_schwarz.cpp @@ -218,7 +218,7 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * // Compute relative error after iteration. double error1 = 0.0; - int norm = 0.0; + double norm = 0.0; for (int m = 0; m < sub_solvers[k]->numSub; m++) { double subdomain_error, subdomain_norm; @@ -243,7 +243,7 @@ void SteadyNSSolver::SchwarzROM(const int M, const int N, ParameterizedProblem * } } // for (int k = 0; k < Ns * Ns; k++) - printf("Iteration %d error: %.4e\n", iter+1, error); + printf("\nIteration %d error: %.4e\n\n", iter+1, error); // Exit the iterations if error is below threshold. if (error <= threshold) { From 0fb124163b1eb57deb265e38f5a22363bb1fdc2f Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Mon, 19 Jan 2026 14:15:52 -0800 Subject: [PATCH 27/27] Schwarz initializes mumps multiple times. delete mumps if it exists. --- src/rom_handler.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index 0e5d0417..c1b75625 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -662,6 +662,7 @@ void MFEMROMHandler::NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec Solver *J_solver = NULL; if (linsol_type == SolverType::DIRECT) { + if (mumps) delete mumps; mumps = new MUMPSSolver(MPI_COMM_SELF); mumps->SetMatrixSymType(mat_type); mumps->SetPrintLevel(jac_print_level);