From cb656e86f51f74cd5e5b3e2d512e99bcd1b8948f Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 14 Mar 2024 10:34:28 -0700 Subject: [PATCH 01/14] initial loading of advection-diffusion solver --- CMakeLists.txt | 3 ++ include/advdiff_solver.hpp | 36 +++++++++++++++++++++ include/poisson_solver.hpp | 3 ++ src/advdiff_solver.cpp | 65 ++++++++++++++++++++++++++++++++++++++ src/poisson_solver.cpp | 12 +++++-- 5 files changed, 116 insertions(+), 3 deletions(-) create mode 100644 include/advdiff_solver.hpp create mode 100644 src/advdiff_solver.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index bffe840f..d09d3d66 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -209,6 +209,9 @@ set(scaleupROMObj_SOURCES include/steady_ns_solver.hpp src/steady_ns_solver.cpp + include/advdiff_solver.hpp + src/advdiff_solver.cpp + include/input_parser.hpp src/input_parser.cpp diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp new file mode 100644 index 00000000..305bd453 --- /dev/null +++ b/include/advdiff_solver.hpp @@ -0,0 +1,36 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#ifndef SCALEUPROM_ADVDIFF_SOLVER_HPP +#define SCALEUPROM_ADVDIFF_SOLVER_HPP + +#include "poisson_solver.hpp" +#include "stokes_solver.hpp" + +// By convention we only use mfem namespace as default, not CAROM. +using namespace mfem; + +class AdvDiffSolver : public PoissonSolver +{ + +friend class ParameterizedProblem; + +protected: + double Pe = 0.0; // Peclet number + +public: + AdvDiffSolver(); + + virtual ~AdvDiffSolver(); + + void BuildDomainOperators() override; + + // Component-wise assembly + void BuildCompROMElement(Array &fes_comp) override; + +protected: + void SetMUMPSSolver() override; +}; + +#endif diff --git a/include/poisson_solver.hpp b/include/poisson_solver.hpp index 2213571b..bf61edb4 100644 --- a/include/poisson_solver.hpp +++ b/include/poisson_solver.hpp @@ -99,6 +99,9 @@ friend class ParameterizedProblem; void SanityCheckOnCoeffs(); virtual void SetParameterizedProblem(ParameterizedProblem *problem); + +protected: + virtual void SetMUMPSSolver(); }; #endif diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp new file mode 100644 index 00000000..3e0f0545 --- /dev/null +++ b/src/advdiff_solver.cpp @@ -0,0 +1,65 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include "advdiff_solver.hpp" +#include "input_parser.hpp" +#include "linalg_utils.hpp" +#include "etc.hpp" + +using namespace std; +using namespace mfem; + +AdvDiffSolver::AdvDiffSolver() + : PoissonSolver() +{ + Pe = config.GetOption("adv-diff/peclet_number", 0.0); +} + +AdvDiffSolver::~AdvDiffSolver() +{ +} + +void AdvDiffSolver::BuildDomainOperators() +{ + PoissonSolver::BuildDomainOperators(); + + for (int m = 0; m < numSub; m++) + { + // as[m]->AddDomainIntegrator(new DiffusionIntegrator); + } +} + +void AdvDiffSolver::BuildCompROMElement(Array &fes_comp) +{ + assert(train_mode == UNIVERSAL); + assert(rom_handler->BasisLoaded()); + + const int num_comp = fes_comp.Size(); + assert(comp_mats.Size() == num_comp); + + for (int c = 0; c < num_comp; c++) + { + Mesh *comp = topol_handler->GetComponentMesh(c); + BilinearForm a_comp(fes_comp[c]); + + a_comp.AddDomainIntegrator(new DiffusionIntegrator); + if (full_dg) + a_comp.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(sigma, kappa)); + + a_comp.Assemble(); + a_comp.Finalize(); + + // Poisson equation has only one solution variable. + comp_mats[c]->SetSize(1, 1); + (*comp_mats[c])(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); + } +} + +void AdvDiffSolver::SetMUMPSSolver() +{ + assert(globalMat_hypre); + mumps = new MUMPSSolver(); + mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); + mumps->SetOperator(*globalMat_hypre); +} \ No newline at end of file diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index 297a9a42..f719ed09 100644 --- a/src/poisson_solver.cpp +++ b/src/poisson_solver.cpp @@ -332,9 +332,7 @@ void PoissonSolver::AssembleOperator() sys_row_starts[1] = globalMat_mono->NumRows(); globalMat_hypre = new HypreParMatrix(MPI_COMM_SELF, sys_glob_size, sys_row_starts, globalMat_mono); - mumps = new MUMPSSolver(); - mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE); - mumps->SetOperator(*globalMat_hypre); + if (direct_solve) SetMUMPSSolver(); } } @@ -586,4 +584,12 @@ void PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem) AddRHSFunction(*(problem->scalar_rhs_ptr)); else AddRHSFunction(0.0); +} + +void PoissonSolver::SetMUMPSSolver() +{ + assert(globalMat_hypre); + mumps = new MUMPSSolver(); + mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE); + mumps->SetOperator(*globalMat_hypre); } \ No newline at end of file From be1ab40546fd580a3fa10d327ae09728fcd77bd5 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 14 Mar 2024 16:35:31 -0700 Subject: [PATCH 02/14] AdvDiffSolver verified with direct solver. --- include/advdiff_solver.hpp | 13 +++ src/advdiff_solver.cpp | 34 +++++++- test/CMakeLists.txt | 1 + test/advdiff_dd_mms.cpp | 68 ++++++++++++++++ test/inputs/dd_mms.yml | 10 ++- test/mms_suite.cpp | 161 +++++++++++++++++++++++++++++++++++++ test/mms_suite.hpp | 25 ++++++ 7 files changed, 309 insertions(+), 3 deletions(-) create mode 100644 test/advdiff_dd_mms.cpp diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index 305bd453..e5cfa7e2 100644 --- a/include/advdiff_solver.hpp +++ b/include/advdiff_solver.hpp @@ -19,6 +19,15 @@ friend class ParameterizedProblem; protected: double Pe = 0.0; // Peclet number + // coefficient for prescribed velocity field. + // can be analytic function or solution from Stokes/SteadyNS equation. + Array flow_coeffs; + + /* + flow solver to obtain the prescribed velocity field. both StokesSolver / SteadyNSSolver can be used. + */ + StokesSolver *stokes_solver = NULL; + public: AdvDiffSolver(); @@ -29,6 +38,10 @@ friend class ParameterizedProblem; // Component-wise assembly void BuildCompROMElement(Array &fes_comp) override; + void SetFlowAtSubdomain(std::function F, const int m=-1); + + void SetParameterizedProblem(ParameterizedProblem *problem) override; + protected: void SetMUMPSSolver() override; }; diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 3e0f0545..03e7f54e 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -13,11 +13,17 @@ using namespace mfem; AdvDiffSolver::AdvDiffSolver() : PoissonSolver() { + // ConvectionIntegrator does not support L2 space. + assert(!full_dg); + Pe = config.GetOption("adv-diff/peclet_number", 0.0); + + flow_coeffs.SetSize(numSub); } AdvDiffSolver::~AdvDiffSolver() { + DeletePointers(flow_coeffs); } void AdvDiffSolver::BuildDomainOperators() @@ -26,12 +32,15 @@ void AdvDiffSolver::BuildDomainOperators() for (int m = 0; m < numSub; m++) { - // as[m]->AddDomainIntegrator(new DiffusionIntegrator); + assert(flow_coeffs[m]); + as[m]->AddDomainIntegrator(new ConvectionIntegrator(*flow_coeffs[m], Pe)); } } void AdvDiffSolver::BuildCompROMElement(Array &fes_comp) { + mfem_error("AdvDiffSolver::BuildCompROMElement is not implemented yet!\n"); + assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); @@ -47,6 +56,8 @@ void AdvDiffSolver::BuildCompROMElement(Array &fes_comp) if (full_dg) a_comp.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(sigma, kappa)); + // TODO(kevin): ConvectionIntegrator needs to be added + a_comp.Assemble(); a_comp.Finalize(); @@ -56,6 +67,27 @@ void AdvDiffSolver::BuildCompROMElement(Array &fes_comp) } } +void AdvDiffSolver::SetFlowAtSubdomain(std::function F, const int m) +{ + assert(flow_coeffs.Size() == numSub); + assert((m == -1) || ((m >= 0) && (m < numSub))); + + if (m == -1) + { + for (int k = 0; k < numSub; k++) + flow_coeffs[k] = new VectorFunctionCoefficient(dim, F); + } + else + flow_coeffs[m] = new VectorFunctionCoefficient(dim, F); +} + +void AdvDiffSolver::SetParameterizedProblem(ParameterizedProblem *problem) +{ + PoissonSolver::SetParameterizedProblem(problem); + + // TODO(kevin): add flow field setup. +} + void AdvDiffSolver::SetMUMPSSolver() { assert(globalMat_hypre); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6c2f3775..9104ee40 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -14,6 +14,7 @@ add_executable(poisson_dd_mms poisson_dd_mms.cpp $ add_executable(stokes_dd_mms stokes_dd_mms.cpp $ $) add_executable(steady_ns_dd_mms steady_ns_dd_mms.cpp $ $) add_executable(linelast_dd_mms linelast_dd_mms.cpp $ $) +add_executable(advdiff_dd_mms advdiff_dd_mms.cpp $ $) add_executable(dg_integ_mms dg_integ_mms.cpp $ $) file(COPY inputs/dd_mms.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs) file(COPY meshes/dd_mms.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes) diff --git a/test/advdiff_dd_mms.cpp b/test/advdiff_dd_mms.cpp new file mode 100644 index 00000000..9dc0b32e --- /dev/null +++ b/test/advdiff_dd_mms.cpp @@ -0,0 +1,68 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include +#include "mms_suite.hpp" + +using namespace std; +using namespace mfem; +using namespace mms::advdiff; + +/** + * Simple smoke test to make sure Google Test is properly linked + */ +TEST(GoogleTestFramework, GoogleTestFrameworkFound) { + SUCCEED(); +} + +// TEST(DDSerialTest, Test_convergence) +// { +// config = InputParser("inputs/dd_mms.yml"); +// CheckConvergence(); + +// return; +// } + +TEST(DDSerialTest, Test_direct_solver) +{ + config = InputParser("inputs/dd_mms.yml"); + config.dict_["solver"]["direct_solve"] = true; + CheckConvergence(); + + return; +} + +// TEST(DDSerial_component_wise_test, Test_convergence) +// { +// config = InputParser("inputs/dd_mms.component.yml"); +// CheckConvergence(); + +// return; +// } + +// TEST(DDSerial_component_3D_hex_test, Test_convergence) +// { +// config = InputParser("inputs/dd_mms.comp.3d.yml"); +// CheckConvergence(); + +// return; +// } + +// TEST(DDSerial_component_3D_tet_test, Test_convergence) +// { +// config = InputParser("inputs/dd_mms.comp.3d.yml"); +// config.dict_["mesh"]["component-wise"]["components"][0]["file"] = "meshes/dd_mms.3d.tet.mesh"; +// CheckConvergence(); + +// return; +// } + +int main(int argc, char* argv[]) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; +} \ No newline at end of file diff --git a/test/inputs/dd_mms.yml b/test/inputs/dd_mms.yml index fb255ad6..6d68543a 100644 --- a/test/inputs/dd_mms.yml +++ b/test/inputs/dd_mms.yml @@ -9,8 +9,11 @@ discretization: full-discrete-galerkin: false visualization: - enable: false - output_dir: dd_mms_output + enabled: false + file_path: + prefix: dd_mms_output + visualize_error: true + unified_paraview: true manufactured_solution: number_of_refinement: 4 @@ -18,5 +21,8 @@ manufactured_solution: amp2: 3.13 constant: 1.11 +adv-diff: + peclet_number: 1.1 + stokes: nu: 1.2 diff --git a/test/mms_suite.cpp b/test/mms_suite.cpp index 2d2ee418..5d95222f 100644 --- a/test/mms_suite.cpp +++ b/test/mms_suite.cpp @@ -597,6 +597,167 @@ namespace linelast } // namespace linelast +namespace advdiff +{ + +double ExactSolution(const Vector &x) +{ + double result = constant; + for (int d = 0; d < x.Size(); d++) + result += amp[d] * sin(2.0 * pi / L[d] * (x(d) - offset[d])); + return result; +} + +void ExactFlow(const Vector &x, Vector &y) +{ + y.SetSize(x.Size()); + y = 0.0; + y(0) = u0; y(1) = v0; + + double xi = 2.0 * pi * wn * (x(0) - uoffset[0]); + double yi = 2.0 * pi * wn * (x(1) - uoffset[1]); + + // incompressible flow field + y(0) += du * cos(xi) * sin(yi); + y(1) += -du * sin(xi) * cos(yi); + + return; +} + +double ExactRHS(const Vector &x) +{ + Vector flow; + ExactFlow(x, flow); + + double result = 0.0; + for (int d = 0; d < x.Size(); d++) + { + result += amp[d] * (2.0 * pi / L[d]) * (2.0 * pi / L[d]) * sin(2.0 * pi / L[d] * (x(d) - offset[d])); + result += Pe * flow(d) * amp[d] * (2.0 * pi / L[d]) * cos(2.0 * pi / L[d] * (x(d) - offset[d])); + } + return result; +} + +AdvDiffSolver *SolveWithRefinement(const int num_refinement) +{ + config.dict_["mesh"]["uniform_refinement"] = num_refinement; + AdvDiffSolver *test = new AdvDiffSolver(); + + test->InitVariables(); + test->InitVisualization(); + + test->AddBCFunction(ExactSolution); + test->SetBdrType(BoundaryType::DIRICHLET); + test->AddRHSFunction(ExactRHS); + test->SetFlowAtSubdomain(ExactFlow); + + FunctionCoefficient exact_sol(ExactSolution); + for (int k = 0; k < test->GetNumSubdomains(); k++) + { + GridFunction *uk = test->GetGridFunction(k); + uk->ProjectCoefficient(exact_sol); + } + BlockVector *exact_U = test->GetSolutionCopy(); + Vector error; + + test->BuildOperators(); + + test->SetupBCOperators(); + + test->Assemble(); + + test->Solve(); + // test->CompareSolution(*exact_U, error); + test->SaveVisualization(); + + delete exact_U; + return test; +} + +void CheckConvergence() +{ + amp[0] = config.GetOption("manufactured_solution/amp1", 0.22); + amp[1] = config.GetOption("manufactured_solution/amp2", 0.13); + amp[2] = config.GetOption("manufactured_solution/amp3", 0.37); + L[0] = config.GetOption("manufactured_solution/L1", 0.31); + L[1] = config.GetOption("manufactured_solution/L2", 0.72); + L[2] = config.GetOption("manufactured_solution/L2", 0.47); + offset[0] = config.GetOption("manufactured_solution/offset1", 0.35); + offset[1] = config.GetOption("manufactured_solution/offset2", 0.73); + offset[2] = config.GetOption("manufactured_solution/offset3", 0.59); + constant = config.GetOption("manufactured_solution/constant", -0.27); + + Pe = config.GetOption("adv-diff/peclet_number", 0.2); + u0 = config.GetOption("manufactured_solution/u0", 1.2); + v0 = config.GetOption("manufactured_solution/v0", -0.7); + du = config.GetOption("manufactured_solution/du", 0.21); + wn = config.GetOption("manufactured_solution/wn", 0.8); + uoffset[0] = config.GetOption("manufactured_solution/uoffset1", 0.51); + uoffset[1] = config.GetOption("manufactured_solution/uoffset2", 0.19); + uoffset[2] = config.GetOption("manufactured_solution/uoffset3", 0.91); + + int num_refine = config.GetOption("manufactured_solution/number_of_refinement", 3); + int base_refine = config.GetOption("manufactured_solution/baseline_refinement", 0); + + // Compare with exact solution + FunctionCoefficient exact_sol(ExactSolution); + + printf("Num. Elem.\tRelative Error\tConvergence Rate\tNorm\n"); + + Vector conv_rate(num_refine); + conv_rate = 0.0; + double error1 = 0.0; + for (int r = base_refine; r < num_refine; r++) + { + AdvDiffSolver *test = SolveWithRefinement(r); + + int order = test->GetDiscretizationOrder(); + int order_quad = max(2, 2*order+1); + const IntegrationRule *irs[Geometry::NumGeom]; + for (int i=0; i < Geometry::NumGeom; ++i) + { + irs[i] = &(IntRules.Get(i, order_quad)); + } + + int numEl = 0; + double norm = 0.0; + for (int k = 0; k < test->GetNumSubdomains(); k++) + { + Mesh *mk = test->GetMesh(k); + norm += pow(ComputeLpNorm(2.0, exact_sol, *mk, irs), 2); + numEl += mk->GetNE(); + } + norm = sqrt(norm); + + double error = 0.0; + for (int k = 0; k < test->GetNumSubdomains(); k++) + { + GridFunction *uk = test->GetGridFunction(k); + error += pow(uk->ComputeLpError(2, exact_sol), 2); + } + error = sqrt(error); + error /= norm; + + if (r > base_refine) + { + conv_rate(r) = error1 / error; + } + printf("%d\t%.15E\t%.15E\t%.15E\n", numEl, error, conv_rate(r), norm); + + // reported convergence rate + if (r > base_refine) + EXPECT_TRUE(conv_rate(r) > pow(2.0, order+1) - 0.5); + + error1 = error; + + delete test; + } + + return; +} + +} // namespace advdiff + namespace fem { diff --git a/test/mms_suite.hpp b/test/mms_suite.hpp index 01b848f0..0e49fc23 100644 --- a/test/mms_suite.hpp +++ b/test/mms_suite.hpp @@ -10,6 +10,7 @@ #include "stokes_solver.hpp" #include "steady_ns_solver.hpp" #include "linelast_solver.hpp" +#include "advdiff_solver.hpp" #include #include #include @@ -90,6 +91,30 @@ void CheckConvergence(); } // namespace linelast +namespace advdiff +{ + +static const double pi = 4.0 * atan(1.0); +// solution parameters +static double amp[3]; +static double L[3]; +static double offset[3]; +static double constant; +// flow parameters +static double Pe; +static double du; +static double wn; +static double uoffset[3]; +static double u0, v0; + +double ExactSolution(const Vector &); +void ExactFlow(const Vector &, Vector &); +double ExactRHS(const Vector &); +AdvDiffSolver *SolveWithRefinement(const int num_refinement); +void CheckConvergence(); + +} // namespace advdiff + namespace fem { From a5bf18ce1340f97b3007d0e1f4f41097d677051f Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 14 Mar 2024 16:47:38 -0700 Subject: [PATCH 03/14] advection diffusion verified with mms. --- include/advdiff_solver.hpp | 2 + src/advdiff_solver.cpp | 77 ++++++++++++++++++++++++++++++++ test/advdiff_dd_mms.cpp | 36 +++++++-------- test/inputs/dd_mms.comp.3d.yml | 6 +++ test/inputs/dd_mms.component.yml | 5 ++- 5 files changed, 107 insertions(+), 19 deletions(-) diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index e5cfa7e2..8d0b18ac 100644 --- a/include/advdiff_solver.hpp +++ b/include/advdiff_solver.hpp @@ -38,6 +38,8 @@ friend class ParameterizedProblem; // Component-wise assembly void BuildCompROMElement(Array &fes_comp) override; + bool Solve() override; + void SetFlowAtSubdomain(std::function F, const int m=-1); void SetParameterizedProblem(ParameterizedProblem *problem) override; diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 03e7f54e..07b83c38 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -67,6 +67,83 @@ void AdvDiffSolver::BuildCompROMElement(Array &fes_comp) } } +bool AdvDiffSolver::Solve() +{ + // If using direct solver, returns always true. + bool converged = true; + + 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); + + // TODO: need to change when the actual parallelization is implemented. + if (direct_solve) + { + assert(mumps); + mumps->SetPrintLevel(print_level); + mumps->Mult(*RHS, *U); + } + else + { + GMRESSolver *solver = NULL; + HypreBoomerAMG *M = NULL; + BlockDiagonalPreconditioner *globalPrec = NULL; + + // HypreBoomerAMG makes a meaningful difference in computation time. + if (use_amg) + { + // Initializating HypreParMatrix needs the monolithic sparse matrix. + assert(globalMat_mono != NULL); + + solver = new GMRESSolver(MPI_COMM_SELF); + + M = new HypreBoomerAMG(*globalMat_hypre); + M->SetPrintLevel(print_level); + + solver->SetPreconditioner(*M); + solver->SetOperator(*globalMat_hypre); + } + else + { + solver = new GMRESSolver(); + + if (config.GetOption("solver/block_diagonal_preconditioner", true)) + { + globalPrec = new BlockDiagonalPreconditioner(var_offsets); + solver->SetPreconditioner(*globalPrec); + } + solver->SetOperator(*globalMat); + } + solver->SetAbsTol(atol); + solver->SetRelTol(rtol); + solver->SetMaxIter(maxIter); + solver->SetPrintLevel(print_level); + + *U = 0.0; + // The time for the setup above is much smaller than this Mult(). + // StopWatch test; + // test.Start(); + solver->Mult(*RHS, *U); + // test.Stop(); + // printf("test: %f seconds.\n", test.RealTime()); + converged = solver->GetConverged(); + + // delete the created objects. + if (use_amg) + { + delete M; + } + else + { + if (globalPrec != NULL) delete globalPrec; + } + delete solver; + } + + return converged; +} + void AdvDiffSolver::SetFlowAtSubdomain(std::function F, const int m) { assert(flow_coeffs.Size() == numSub); diff --git a/test/advdiff_dd_mms.cpp b/test/advdiff_dd_mms.cpp index 9dc0b32e..9914de13 100644 --- a/test/advdiff_dd_mms.cpp +++ b/test/advdiff_dd_mms.cpp @@ -16,13 +16,13 @@ TEST(GoogleTestFramework, GoogleTestFrameworkFound) { SUCCEED(); } -// TEST(DDSerialTest, Test_convergence) -// { -// config = InputParser("inputs/dd_mms.yml"); -// CheckConvergence(); +TEST(DDSerialTest, Test_convergence) +{ + config = InputParser("inputs/dd_mms.yml"); + CheckConvergence(); -// return; -// } + return; +} TEST(DDSerialTest, Test_direct_solver) { @@ -33,21 +33,21 @@ TEST(DDSerialTest, Test_direct_solver) return; } -// TEST(DDSerial_component_wise_test, Test_convergence) -// { -// config = InputParser("inputs/dd_mms.component.yml"); -// CheckConvergence(); +TEST(DDSerial_component_wise_test, Test_convergence) +{ + config = InputParser("inputs/dd_mms.component.yml"); + CheckConvergence(); -// return; -// } + return; +} -// TEST(DDSerial_component_3D_hex_test, Test_convergence) -// { -// config = InputParser("inputs/dd_mms.comp.3d.yml"); -// CheckConvergence(); +TEST(DDSerial_component_3D_hex_test, Test_convergence) +{ + config = InputParser("inputs/dd_mms.comp.3d.yml"); + CheckConvergence(); -// return; -// } + return; +} // TEST(DDSerial_component_3D_tet_test, Test_convergence) // { diff --git a/test/inputs/dd_mms.comp.3d.yml b/test/inputs/dd_mms.comp.3d.yml index 44b668f7..529e4c55 100644 --- a/test/inputs/dd_mms.comp.3d.yml +++ b/test/inputs/dd_mms.comp.3d.yml @@ -39,6 +39,9 @@ discretization: order: 2 full-discrete-galerkin: false +solver: + direct_solve: true + visualization: enable: false output_dir: dd_mms_output @@ -50,5 +53,8 @@ manufactured_solution: amp3: 2.77 constant: 1.11 +adv-diff: + peclet_number: 1.1 + stokes: nu: 1.2 diff --git a/test/inputs/dd_mms.component.yml b/test/inputs/dd_mms.component.yml index 25d4421b..69c1797d 100644 --- a/test/inputs/dd_mms.component.yml +++ b/test/inputs/dd_mms.component.yml @@ -39,5 +39,8 @@ manufactured_solution: amp2: 3.13 constant: 1.11 +adv-diff: + peclet_number: 1.1 + stokes: - nu: 1.2 \ No newline at end of file + nu: 1.2 From a3eccfefa27cd62a5e4a38b35bfe15c3943dd614 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 14 Mar 2024 16:52:43 -0700 Subject: [PATCH 04/14] component-wise verification with mms for advection-diffusion. --- test/gmsh/multi_comp_dd_mms.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/gmsh/multi_comp_dd_mms.cpp b/test/gmsh/multi_comp_dd_mms.cpp index f956c485..a755c113 100644 --- a/test/gmsh/multi_comp_dd_mms.cpp +++ b/test/gmsh/multi_comp_dd_mms.cpp @@ -16,7 +16,7 @@ TEST(GoogleTestFramework, GoogleTestFrameworkFound) { SUCCEED(); } -TEST(DDSerialTest, Test_convergence) +TEST(Poisson, Test_convergence) { config = InputParser("test.component.yml"); CheckConvergence(); @@ -24,6 +24,15 @@ TEST(DDSerialTest, Test_convergence) return; } +TEST(AdvDiff, Test_convergence) +{ + config = InputParser("test.component.yml"); + config.dict_["adv-diff"]["peclet_number"] = 1.1; + mms::advdiff::CheckConvergence(); + + return; +} + int main(int argc, char* argv[]) { MPI_Init(&argc, &argv); From 1c08528f723eb8a1ef3619e3a8a42c72bb7c1c0b Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 15 Mar 2024 13:12:13 -0700 Subject: [PATCH 05/14] AdvDiffSolver::GetFlowField and AdvDiffFlowPastArray class. still needs debugging. --- include/advdiff_solver.hpp | 6 +++ include/parameterized_problem.hpp | 68 +++++++++++++++++++++++--- src/advdiff_solver.cpp | 47 +++++++++++++++++- src/main_workflow.cpp | 2 + src/parameterized_problem.cpp | 78 +++++++++++++++++++++++++++++ src/poisson_solver.cpp | 3 ++ src/stokes_solver.cpp | 3 ++ test/CMakeLists.txt | 1 + test/inputs/advdiff.base.yml | 81 +++++++++++++++++++++++++++++++ test/test_workflow.cpp | 30 ++++++++++++ 10 files changed, 310 insertions(+), 9 deletions(-) create mode 100644 test/inputs/advdiff.base.yml diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index 8d0b18ac..303e7662 100644 --- a/include/advdiff_solver.hpp +++ b/include/advdiff_solver.hpp @@ -27,6 +27,9 @@ friend class ParameterizedProblem; flow solver to obtain the prescribed velocity field. both StokesSolver / SteadyNSSolver can be used. */ StokesSolver *stokes_solver = NULL; + bool load_flow = false; + bool save_flow = false; + std::string flow_file = ""; public: AdvDiffSolver(); @@ -46,6 +49,9 @@ friend class ParameterizedProblem; protected: void SetMUMPSSolver() override; + +private: + void GetFlowField(ParameterizedProblem *flow_problem); }; #endif diff --git a/include/parameterized_problem.hpp b/include/parameterized_problem.hpp index e815a3e2..dccca6f0 100644 --- a/include/parameterized_problem.hpp +++ b/include/parameterized_problem.hpp @@ -92,9 +92,24 @@ extern double rforce_f; void tip_force(const Vector &x, Vector &u); } +namespace advdiff_problem +{ + +extern bool analytic_flow; + +namespace advdiff_flow_past_array +{ + extern double q0, dq, qoffset; + extern Vector qk; + + double qbdr(const Vector &x); +} // namespace advdiff_flow_past_array + +} // namespace advdiff_problem + } -enum BoundaryType +enum class BoundaryType { ZERO, DIRICHLET, @@ -169,8 +184,8 @@ class PoissonComponent : public PoissonProblem public: PoissonComponent(); virtual ~PoissonComponent() {}; - virtual void SetParams(const std::string &key, const double &value) override; - virtual void SetParams(const Array &indexes, const Vector &values) override; + void SetParams(const std::string &key, const double &value) override; + void SetParams(const Array &indexes, const Vector &values) override; private: void SetBattr(); @@ -206,15 +221,17 @@ class StokesComponent : public StokesProblem class StokesFlowPastArray : public StokesComponent { +friend class AdvDiffFlowPastArray; + public: StokesFlowPastArray(); - virtual void SetParams(const std::string &key, const double &value); - virtual void SetParams(const Array &indexes, const Vector &values); + virtual void SetParams(const std::string &key, const double &value) override; + virtual void SetParams(const Array &indexes, const Vector &values) override; -private: +protected: Vector *u0; - void SetBattr(); + virtual void SetBattr(); }; class LinElastProblem : public ParameterizedProblem @@ -251,6 +268,43 @@ class LinElastForceCantilever : public LinElastProblem LinElastForceCantilever(); }; +namespace function_factory +{ + +namespace advdiff_problem +{ +/* + flow_problem will be passed down to StokesSolver/SteadyNSSolver for obtaining velocity field. + It must be set appropriately within each AdvDiffSolver problems. +*/ +extern StokesProblem *flow_problem; + +} // namespace advdiff_problem + +} // namespace function_factory + +class AdvDiffFlowPastArray : public StokesFlowPastArray +{ +protected: + /* + flow_problem shares the same pointers with this class. + Thus every parameter set by this class is reflected to StokesFlowPastArrayProblem as well. + flow_problem will be passed down to StokesSolver/SteadyNSSolver for obtaining velocity field. + */ + StokesFlowPastArray *flow_problem = NULL; + +public: + AdvDiffFlowPastArray(); + virtual ~AdvDiffFlowPastArray(); + +protected: + void SetBattr() override + { + StokesFlowPastArray::SetBattr(); + flow_problem->SetBattr(); + } +}; + ParameterizedProblem* InitParameterizedProblem(); #endif diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 07b83c38..e100c89a 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -19,11 +19,17 @@ AdvDiffSolver::AdvDiffSolver() Pe = config.GetOption("adv-diff/peclet_number", 0.0); flow_coeffs.SetSize(numSub); + + save_flow = config.GetOption("adv-diff/save_flow", false); + load_flow = config.GetOption("adv-diff/load_flow", false); + if (save_flow || load_flow) + flow_file = config.GetRequiredOption("adv-diff/flow_file"); } AdvDiffSolver::~AdvDiffSolver() { DeletePointers(flow_coeffs); + delete stokes_solver; } void AdvDiffSolver::BuildDomainOperators() @@ -160,9 +166,10 @@ void AdvDiffSolver::SetFlowAtSubdomain(std::functionSetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); mumps->SetOperator(*globalMat_hypre); +} + +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(); + stokes_solver->InitVisualization(); + + if (load_flow && FileExists(flow_file)) + stokes_solver->LoadSolution(flow_file); + else + { + stokes_solver->SetParameterizedProblem(flow_problem); + // currently only support FOM. + stokes_solver->BuildOperators(); + stokes_solver->SetupBCOperators(); + stokes_solver->Assemble(); + stokes_solver->Solve(); + } + + if (save_flow) + stokes_solver->SaveSolution(flow_file); + + DeletePointers(flow_coeffs); + const int stokes_numvar = stokes_solver->GetNumVar(); + for (int m = 0; m < numSub; m++) + flow_coeffs[m] = new VectorGridFunctionCoefficient(stokes_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(). + */ } \ No newline at end of file diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index 56b1a675..08756d63 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -9,6 +9,7 @@ #include "linelast_solver.hpp" #include "stokes_solver.hpp" #include "steady_ns_solver.hpp" +#include "advdiff_solver.hpp" #include "etc.hpp" #include #include @@ -55,6 +56,7 @@ MultiBlockSolver* InitSolver() else if (solver_type == "stokes") { solver = new StokesSolver; } else if (solver_type == "steady-ns") { solver = new SteadyNSSolver; } else if (solver_type == "linelast") { solver = new LinElastSolver; } + else if (solver_type == "adv-diff") { solver = new AdvDiffSolver; } else { printf("Unknown MultiBlockSolver %s!\n", solver_type.c_str()); diff --git a/src/parameterized_problem.cpp b/src/parameterized_problem.cpp index a40119e2..59a7b73c 100644 --- a/src/parameterized_problem.cpp +++ b/src/parameterized_problem.cpp @@ -200,6 +200,32 @@ void tip_force(const Vector &x, Vector &f) } // namespace linelast_force +namespace advdiff_problem +{ + +bool analytic_flow; +StokesProblem *flow_problem; + +namespace advdiff_flow_past_array +{ + +double q0, dq, qoffset; +Vector qk; + +double qbdr(const Vector &x) +{ + assert(qk.Size() >= x.Size()); + double tmp = 0.0; + for (int d = 0; d < x.Size(); d++) + tmp += advdiff_flow_past_array::qk(d) * x(d); + tmp += advdiff_flow_past_array::qoffset; + return advdiff_flow_past_array::q0 + advdiff_flow_past_array::dq * sin(2.0 * pi * tmp); +} + +} // namespace advdiff_flow_past_array + +} // namespace advdiff_problem + } // namespace function_factory ParameterizedProblem::ParameterizedProblem() @@ -303,6 +329,10 @@ ParameterizedProblem* InitParameterizedProblem() { problem = new LinElastForceCantilever(); } + else if (problem_name == "advdiff_flow_past_array") + { + problem = new AdvDiffFlowPastArray(); + } else { mfem_error("Unknown parameterized problem name!\n"); @@ -762,3 +792,51 @@ LinElastForceCantilever::LinElastForceCantilever() general_vector_ptr.SetSize(1); general_vector_ptr[0] = NULL; } + +/* + AdvDiffFlowPastArray +*/ + +AdvDiffFlowPastArray::AdvDiffFlowPastArray() + : StokesFlowPastArray(), flow_problem(new StokesFlowPastArray) +{ + function_factory::advdiff_problem::analytic_flow = false; + function_factory::advdiff_problem::flow_problem = flow_problem; + + bdr_type.SetSize(5); + bdr_type = BoundaryType::DIRICHLET; + bdr_type[4] = BoundaryType::NEUMANN; + + scalar_rhs_ptr = NULL; + scalar_bdr_ptr.SetSize(5); + scalar_bdr_ptr = &(function_factory::advdiff_problem::advdiff_flow_past_array::qbdr); + + // q0 + dq + qoffset + qk(3) + param_num += 1 + 1 + 1 + 3; + const int p0 = flow_problem->GetNumParams(); + + param_map["q0"] = p0; + param_map["dq"] = p0 + 1; + param_map["qoffset"] = p0 + 2; + param_map["qk_x"] = p0 + 3; + param_map["qk_y"] = p0 + 4; + param_map["qk_z"] = p0 + 5; + + // default values. + function_factory::advdiff_problem::advdiff_flow_past_array::q0 = 1.0; + function_factory::advdiff_problem::advdiff_flow_past_array::dq = 0.1; + function_factory::advdiff_problem::advdiff_flow_past_array::qoffset = 0.0; + function_factory::advdiff_problem::advdiff_flow_past_array::qk.SetSize(3); + function_factory::advdiff_problem::advdiff_flow_past_array::qk = 0.0; + + param_ptr.Append(&(function_factory::advdiff_problem::advdiff_flow_past_array::q0)); + param_ptr.Append(&(function_factory::advdiff_problem::advdiff_flow_past_array::dq)); + param_ptr.Append(&(function_factory::advdiff_problem::advdiff_flow_past_array::qoffset)); + for (int j = 0; j < 3; j++) + param_ptr.Append(&(function_factory::advdiff_problem::advdiff_flow_past_array::qk(j))); +} + +AdvDiffFlowPastArray::~AdvDiffFlowPastArray() +{ + delete flow_problem; +} diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index f719ed09..c768b37f 100644 --- a/src/poisson_solver.cpp +++ b/src/poisson_solver.cpp @@ -554,6 +554,9 @@ void PoissonSolver::SanityCheckOnCoeffs() void PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem) { + /* set up boundary types */ + MultiBlockSolver::SetParameterizedProblem(problem); + // clean up rhs for parametrized problem. if (rhs_coeffs.Size() > 0) { diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index 65b4cb50..b39a2f8a 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -1078,6 +1078,9 @@ void StokesSolver::SanityCheckOnCoeffs() void StokesSolver::SetParameterizedProblem(ParameterizedProblem *problem) { + /* set up boundary types */ + MultiBlockSolver::SetParameterizedProblem(problem); + nu = function_factory::stokes_problem::nu; delete nu_coeff; nu_coeff = new ConstantCoefficient(nu); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 9104ee40..2f2ec3fd 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -33,6 +33,7 @@ add_executable(test_workflow test_workflow.cpp $) file(COPY inputs/test.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs) file(COPY inputs/stokes.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs) file(COPY inputs/steady_ns.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs) +file(COPY inputs/advdiff.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs) file(COPY meshes/test.2x2.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes) file(COPY meshes/test.2x1.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes) diff --git a/test/inputs/advdiff.base.yml b/test/inputs/advdiff.base.yml new file mode 100644 index 00000000..8c905440 --- /dev/null +++ b/test/inputs/advdiff.base.yml @@ -0,0 +1,81 @@ +main: +#mode: run_example/sample_generation/build_rom/single_run + solver: adv-diff + mode: single_run + use_rom: true + +adv-diff: + peclet_number: 1.1 + save_flow: true + flow_file: advdiff.flow_field.h5 + +mesh: + filename: meshes/test.2x2.mesh + uniform_refinement: 1 + +domain-decomposition: + type: interior_penalty + +discretization: + order: 1 + full-discrete-galerkin: false + +visualization: + enabled: false + unified_paraview: true + file_path: + prefix: sample_gen_output + +parameterized_problem: + name: advdiff_flow_past_array + +single_run: + advdiff_flow_past_array: + nu: 2.5 + u0_x: 1.0 + u0_y: -1.0 + du_x: 0.1 + du_y: 0.1 + k_u_x: 0.5 + k_u_y: 0.8 + k_v_x: 0.7 + k_v_y: 0.2 + q0: 1.0 + dq: 0.8 + qoffset: 0.0 + qk_x: 1.5 + qk_y: 0.9 + +sample_generation: + maximum_number_of_snapshots: 400 + component_sampling: false + file_path: + prefix: "advdiff_flow_past_array" + parameters: + - key: single_run/advdiff_flow_past_array/qk_x + type: double + sample_size: 3 + minimum: 1.0 + maximum: 2.0 + +basis: + prefix: "advdiff0" + number_of_basis: 3 + tags: + - name: comp0 + svd: + save_spectrum: true + update_right_sv: false + visualization: + enabled: false + +model_reduction: + # individual/universal + subdomain_training: individual + save_operator: + level: global + prefix: "proj_inv" + compare_solution: + enabled: true + linear_solver_type: direct + linear_system_type: us diff --git a/test/test_workflow.cpp b/test/test_workflow.cpp index 35fa3617..62b92e32 100644 --- a/test/test_workflow.cpp +++ b/test/test_workflow.cpp @@ -709,6 +709,36 @@ TEST(LinElast_Workflow, ComponentWiseTest) return; } +TEST(AdvDiff_Workflow, MFEMIndividualTest) +{ + config = InputParser("inputs/advdiff.base.yml"); + for (int k = 0; k < 4; k++) + config.dict_["basis"]["tags"][k]["name"] = "dom" + std::to_string(k); + + config.dict_["model_reduction"]["visualization"]["enabled"] = true; + config.dict_["model_reduction"]["visualization"]["prefix"] = "basis_paraview"; + + config.dict_["main"]["mode"] = "sample_generation"; + config.dict_["main"]["use_rom"] = false; + GenerateSamples(MPI_COMM_WORLD); + config.dict_["main"]["use_rom"] = true; + + config.dict_["main"]["mode"] = "train_rom"; + TrainROM(MPI_COMM_WORLD); + + config.dict_["main"]["mode"] = "build_rom"; + BuildROM(MPI_COMM_WORLD); + + config.dict_["main"]["mode"] = "single_run"; + double error = SingleRun(MPI_COMM_WORLD); + + // This reproductive case must have a very small error at the level of finite-precision. + printf("Error: %.15E\n", error); + EXPECT_TRUE(error < threshold); + + return; +} + int main(int argc, char *argv[]) { ::testing::InitGoogleTest(&argc, argv); From 02e50e16120ed18f67d1cc7f7cdcfbdb6c4063fc Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 15 Mar 2024 14:09:48 -0700 Subject: [PATCH 06/14] VectorGridFunctionCoefficient needs mfem v4.6. updated docker. --- docker/Dockerfile | 4 +++- sketches/ns_dg_mms.cpp | 2 +- sketches/ns_rom.cpp | 2 +- src/advdiff_solver.cpp | 3 ++- src/linelast_solver.cpp | 2 +- src/poisson_solver.cpp | 4 ++-- src/rom_handler.cpp | 4 ++-- src/steady_ns_solver.cpp | 2 +- src/stokes_solver.cpp | 4 ++-- 9 files changed, 15 insertions(+), 12 deletions(-) diff --git a/docker/Dockerfile b/docker/Dockerfile index 2272c86e..ba97fd92 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -29,6 +29,8 @@ WORKDIR /usr/lib RUN arch=$(uname -m) && sudo ln -s /usr/lib/${arch}-linux-gnu/libscalapack-openmpi.so.2.1.0 ./ # re-install mfem with mumps, with cmake +WORKDIR $LIB_DIR/mfem +RUN sudo git checkout v4.6 WORKDIR $LIB_DIR/mfem/build RUN sudo -E cmake .. -DBUILD_SHARED_LIBS=YES -DMFEM_USE_MPI=YES -DMFEM_USE_GSLIB=${MG} -DMFEM_USE_METIS=YES -DMFEM_USE_METIS_5=YES -DMFEM_USE_MUMPS=YES -DGSLIB_DIR="$LIB_DIR/gslib/build" -DMUMPS_DIR="$LIB_DIR/mumps/build/local" RUN sudo -E make -j 16 @@ -54,7 +56,7 @@ WORKDIR ./libROM RUN sudo git pull && sudo git checkout mpi-io WORKDIR ./build # libROM is using the MFEM without MUMPS right now. -RUN sudo cmake .. -DCMAKE_TOOLCHAIN_FILE=${TOOLCHAIN_FILE} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DUSE_MFEM=${USE_MFEM} -DMFEM_USE_GSLIB=${MFEM_USE_GSLIB} +RUN sudo cmake .. -DCMAKE_TOOLCHAIN_FILE=${TOOLCHAIN_FILE} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DUSE_MFEM=OFF -DMFEM_USE_GSLIB=${MFEM_USE_GSLIB} RUN sudo make -j 16 ENV LIBROM_DIR=$LIB_DIR/libROM diff --git a/sketches/ns_dg_mms.cpp b/sketches/ns_dg_mms.cpp index e771ee31..45dd2be3 100644 --- a/sketches/ns_dg_mms.cpp +++ b/sketches/ns_dg_mms.cpp @@ -413,7 +413,7 @@ int main(int argc, char *argv[]) MUMPSSolver *J_mumps = NULL; if (direct_solve) { - J_mumps = new MUMPSSolver; + J_mumps = new MUMPSSolver(MPI_COMM_WORLD); J_mumps->SetPrintLevel(-1); J_mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); J_solver = J_mumps; diff --git a/sketches/ns_rom.cpp b/sketches/ns_rom.cpp index e578d030..38f0a110 100644 --- a/sketches/ns_rom.cpp +++ b/sketches/ns_rom.cpp @@ -242,7 +242,7 @@ class SteadyNavierStokes : public Operator if (direct_solve) { - J_mumps = new MUMPSSolver(); + J_mumps = new MUMPSSolver(MPI_COMM_WORLD); J_mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); J_solver = J_mumps; } diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index e100c89a..61bd190d 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -19,6 +19,7 @@ AdvDiffSolver::AdvDiffSolver() Pe = config.GetOption("adv-diff/peclet_number", 0.0); flow_coeffs.SetSize(numSub); + flow_coeffs = NULL; save_flow = config.GetOption("adv-diff/save_flow", false); load_flow = config.GetOption("adv-diff/load_flow", false); @@ -175,7 +176,7 @@ void AdvDiffSolver::SetParameterizedProblem(ParameterizedProblem *problem) void AdvDiffSolver::SetMUMPSSolver() { assert(globalMat_hypre); - mumps = new MUMPSSolver(); + mumps = new MUMPSSolver(MPI_COMM_SELF); mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); mumps->SetOperator(*globalMat_hypre); } diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index d7e7259c..a3e88e58 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -292,7 +292,7 @@ void LinElastSolver::AssembleOperator() sys_row_starts[1] = globalMat_mono->NumRows(); globalMat_hypre = new HypreParMatrix(MPI_COMM_WORLD, sys_glob_size, sys_row_starts, globalMat_mono); - mumps = new MUMPSSolver(); + mumps = new MUMPSSolver(MPI_COMM_SELF); mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE); mumps->SetOperator(*globalMat_hypre); } diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index c768b37f..d9c8844d 100644 --- a/src/poisson_solver.cpp +++ b/src/poisson_solver.cpp @@ -556,7 +556,7 @@ void PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem) { /* set up boundary types */ MultiBlockSolver::SetParameterizedProblem(problem); - + // clean up rhs for parametrized problem. if (rhs_coeffs.Size() > 0) { @@ -592,7 +592,7 @@ void PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem) void PoissonSolver::SetMUMPSSolver() { assert(globalMat_hypre); - mumps = new MUMPSSolver(); + mumps = new MUMPSSolver(MPI_COMM_SELF); mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE); mumps->SetOperator(*globalMat_hypre); } \ No newline at end of file diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index 6f001d4e..28690f87 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -674,7 +674,7 @@ void MFEMROMHandler::NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec Solver *J_solver = NULL; if (linsol_type == SolverType::DIRECT) { - mumps = new MUMPSSolver(); + mumps = new MUMPSSolver(MPI_COMM_SELF); mumps->SetMatrixSymType(mat_type); mumps->SetPrintLevel(jac_print_level); J_solver = mumps; @@ -1086,7 +1086,7 @@ void MFEMROMHandler::SetupDirectSolver() sys_row_starts[1] = romMat_mono->NumRows(); romMat_hypre = new HypreParMatrix(MPI_COMM_SELF, sys_glob_size, sys_row_starts, romMat_mono); - mumps = new MUMPSSolver(); + mumps = new MUMPSSolver(MPI_COMM_SELF); mumps->SetMatrixSymType(mat_type); mumps->SetOperator(*romMat_hypre); } diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index 7282fab0..d75785ab 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -645,7 +645,7 @@ bool SteadyNSSolver::Solve() if (direct_solve) { - mumps = new MUMPSSolver(); + mumps = new MUMPSSolver(MPI_COMM_SELF); mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); mumps->SetPrintLevel(jac_print_level); J_solver = mumps; diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index b39a2f8a..19e0db9d 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -492,7 +492,7 @@ void StokesSolver::SetupMUMPSSolver() sys_row_starts[1] = systemOp_mono->NumRows(); systemOp_hypre = new HypreParMatrix(MPI_COMM_SELF, sys_glob_size, sys_row_starts, systemOp_mono); - mumps = new MUMPSSolver(); + mumps = new MUMPSSolver(MPI_COMM_SELF); mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE); mumps->SetOperator(*systemOp_hypre); } @@ -1080,7 +1080,7 @@ void StokesSolver::SetParameterizedProblem(ParameterizedProblem *problem) { /* set up boundary types */ MultiBlockSolver::SetParameterizedProblem(problem); - + nu = function_factory::stokes_problem::nu; delete nu_coeff; nu_coeff = new ConstantCoefficient(nu); From 1954a010b6d77bc30caf71ce63f7e8b96abc7471 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 15 Mar 2024 16:24:07 -0700 Subject: [PATCH 07/14] reflecting librom update --- src/rom_handler.cpp | 2 +- src/sample_generator.cpp | 7 +++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index 28690f87..068ecd32 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -269,7 +269,7 @@ void ROMHandlerBase::LoadReducedBasis() */ { int local_dim = CAROM::split_dimension(dim_ref_basis[k], MPI_COMM_WORLD); - basis_reader = new CAROM::BasisReader(basis_name + basis_tags[k], CAROM::Database::HDF5_MPIO, local_dim); + basis_reader = new CAROM::BasisReader(basis_name + basis_tags[k], CAROM::Database::formats::HDF5_MPIO, local_dim); carom_ref_basis[k] = new CAROM::Matrix(*basis_reader->getSpatialBasis(num_ref_basis[k])); carom_ref_basis[k]->gather(); diff --git a/src/sample_generator.cpp b/src/sample_generator.cpp index 3ffb3af1..9e99682c 100644 --- a/src/sample_generator.cpp +++ b/src/sample_generator.cpp @@ -186,7 +186,7 @@ void SampleGenerator::AddSnapshotGenerator(const int &fom_vdofs, const std::stri snapshot_options.Append(new CAROM::Options(fom_vdofs, max_num_snapshots, 1, update_right_SV)); snapshot_options.Last()->static_svd_preserve_snapshot = true; - snapshot_generators.Append(new CAROM::BasisGenerator(*(snapshot_options.Last()), incremental, filename, CAROM::Database::HDF5_MPIO)); + snapshot_generators.Append(new CAROM::BasisGenerator(*(snapshot_options.Last()), incremental, filename, CAROM::Database::formats::HDF5_MPIO)); basis_tag2idx[basis_tag] = basis_tags.size(); basis_tags.push_back(basis_tag); @@ -264,7 +264,7 @@ void SampleGenerator::FormReducedBasis(const std::string &basis_prefix, CAROM::BasisGenerator *basis_generator = snapshot_generators.Last(); for (int s = 0; s < file_list.size(); s++) - basis_generator->loadSamples(file_list[s], "snapshot", 1e9, CAROM::Database::HDF5_MPIO); + basis_generator->loadSamples(file_list[s], "snapshot", 1e9, CAROM::Database::formats::HDF5_MPIO); basis_generator->endSamples(); // save the merged basis file SaveSV(basis_generator, basis_name, ref_num_basis); @@ -289,8 +289,7 @@ const int SampleGenerator::GetDimFromSnapshots(const std::string &filename) file_id = H5Fopen(filename_ext.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); assert(file_id >= 0); - // TODO(kevin): will have to reflect dataset name update. - hdf5_utils::ReadDataset(file_id, "snapshot_matrix_num_rows_000000", nrows); + hdf5_utils::ReadDataset(file_id, "snapshot_matrix_num_rows", nrows); assert(nrows[0] > 0); errf = H5Fclose(file_id); From 05ce38bfd73a32e65f76dc9fe851040659618635 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 15 Mar 2024 17:36:43 -0700 Subject: [PATCH 08/14] bug fix for SteadyNSSolver: SetMUMPSSolver does not set operator in AssembleOperator. --- include/steady_ns_solver.hpp | 33 ++++++++++++++++----------------- include/stokes_solver.hpp | 2 +- src/steady_ns_solver.cpp | 19 +++++++------------ src/stokes_solver.cpp | 6 +++--- 4 files changed, 27 insertions(+), 33 deletions(-) diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 595e789c..04009ad8 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -150,35 +150,34 @@ friend class SteadyNSOperator; using StokesSolver::GetVariableNames; - virtual void InitVariables(); + void InitVariables() override; - virtual void BuildOperators() override; - virtual void BuildDomainOperators(); + void BuildDomainOperators() override; - virtual void SetupRHSBCOperators() override; - virtual void SetupDomainBCOperators() override; + void SetupRHSBCOperators() override; + void SetupDomainBCOperators() override; - virtual void Assemble(); + void AssembleOperator() override; - virtual void SaveROMOperator(const std::string input_prefix=""); - virtual void LoadROMOperatorFromFile(const std::string input_prefix=""); + void SaveROMOperator(const std::string input_prefix="") override; + void LoadROMOperatorFromFile(const std::string input_prefix="") override; // Component-wise assembly - virtual void BuildCompROMElement(Array &fes_comp); + void BuildCompROMElement(Array &fes_comp) override; // virtual void BuildBdrROMElement(Array &fes_comp); // virtual void BuildInterfaceROMElement(Array &fes_comp); - virtual void SaveCompBdrROMElement(hid_t &file_id) override; - virtual void LoadCompBdrROMElement(hid_t &file_id) override; + void SaveCompBdrROMElement(hid_t &file_id) override; + void LoadCompBdrROMElement(hid_t &file_id) override; - virtual bool Solve(); + bool Solve() override; - virtual void ProjectOperatorOnReducedBasis(); + void ProjectOperatorOnReducedBasis() override; - virtual void SolveROM() override; + void SolveROM() override; - virtual void TrainEQP(SampleGenerator *sample_generator) override; - virtual void SaveEQP() override; - virtual void LoadEQP() override; + void TrainEQP(SampleGenerator *sample_generator) override; + void SaveEQP() override; + void LoadEQP() override; private: DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace); diff --git a/include/stokes_solver.hpp b/include/stokes_solver.hpp index c1919196..4a577f26 100644 --- a/include/stokes_solver.hpp +++ b/include/stokes_solver.hpp @@ -179,7 +179,7 @@ friend class ParameterizedProblem; // For bilinear case. // system-specific. virtual void AssembleInterfaceMatrices(); - virtual void SetupMUMPSSolver(); + virtual void SetupMUMPSSolver(bool set_oper); virtual void SetupPressureMassMatrix(); // Component-wise assembly diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index d75785ab..ee89d7ba 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -329,13 +329,6 @@ void SteadyNSSolver::InitVariables() } } -void SteadyNSSolver::BuildOperators() -{ - BuildRHSOperators(); - - BuildDomainOperators(); -} - void SteadyNSSolver::BuildDomainOperators() { StokesSolver::BuildDomainOperators(); @@ -427,13 +420,15 @@ void SteadyNSSolver::SetupDomainBCOperators() } -void SteadyNSSolver::Assemble() +void SteadyNSSolver::AssembleOperator() { - StokesSolver::AssembleRHS(); + StokesSolver::AssembleOperatorBase(); - StokesSolver::AssembleOperator(); - - // nonlinear operator? + if (direct_solve) + StokesSolver::SetupMUMPSSolver(false); + else + // pressure mass matrix for preconditioner. + StokesSolver::SetupPressureMassMatrix(); } void SteadyNSSolver::SaveROMOperator(const std::string input_prefix) diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index 19e0db9d..6482c13e 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -399,7 +399,7 @@ void StokesSolver::AssembleOperator() AssembleOperatorBase(); if (direct_solve) - SetupMUMPSSolver(); + SetupMUMPSSolver(true); else // pressure mass matrix for preconditioner. SetupPressureMassMatrix(); @@ -477,7 +477,7 @@ void StokesSolver::AssembleOperatorBase() systemOp->SetBlock(1,0, B); } -void StokesSolver::SetupMUMPSSolver() +void StokesSolver::SetupMUMPSSolver(bool set_oper) { assert(systemOp); delete systemOp_mono; @@ -494,7 +494,7 @@ void StokesSolver::SetupMUMPSSolver() mumps = new MUMPSSolver(MPI_COMM_SELF); mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE); - mumps->SetOperator(*systemOp_hypre); + if (set_oper) mumps->SetOperator(*systemOp_hypre); } void StokesSolver::SetupPressureMassMatrix() From 28a13e883f615caf38bce28337c5ed2d9e8aaadd Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 15 Mar 2024 17:55:18 -0700 Subject: [PATCH 09/14] minor fix in sketches. --- sketches/ns_rom.cpp | 2 +- sketches/stokes_mms4.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/sketches/ns_rom.cpp b/sketches/ns_rom.cpp index 38f0a110..f51a70bf 100644 --- a/sketches/ns_rom.cpp +++ b/sketches/ns_rom.cpp @@ -1768,7 +1768,7 @@ int main(int argc, char *argv[]) J_gmres.SetMaxIter(maxIter); J_gmres.SetPrintLevel(-1); - MUMPSSolver J_mumps; + MUMPSSolver J_mumps(MPI_COMM_WORLD); J_mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); Solver *J_solver; diff --git a/sketches/stokes_mms4.cpp b/sketches/stokes_mms4.cpp index 2bdb114b..0bcc9d28 100644 --- a/sketches/stokes_mms4.cpp +++ b/sketches/stokes_mms4.cpp @@ -322,7 +322,7 @@ int main(int argc, char *argv[]) HYPRE_BigInt sys_row_starts[2] = {0, sys_mat->NumRows()}; HypreParMatrix Aop(MPI_COMM_WORLD, sys_glob_size, sys_row_starts, sys_mat); - MUMPSSolver mumps; + MUMPSSolver mumps(MPI_COMM_WORLD); mumps.SetPrintLevel(1); mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); // mumps.SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE); From 7da525e78229d0bee628f763a0f63a078d56d386 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Fri, 15 Mar 2024 18:14:50 -0700 Subject: [PATCH 10/14] StokesSolver now must use SYMMETRIC_INDEFINITE. --- src/stokes_solver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index 6482c13e..23afb953 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -493,7 +493,7 @@ void StokesSolver::SetupMUMPSSolver(bool set_oper) systemOp_hypre = new HypreParMatrix(MPI_COMM_SELF, sys_glob_size, sys_row_starts, systemOp_mono); mumps = new MUMPSSolver(MPI_COMM_SELF); - mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE); + mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_INDEFINITE); if (set_oper) mumps->SetOperator(*systemOp_hypre); } From 1262a7a8f818f1b48fd464458a202126435d59fb Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Sat, 16 Mar 2024 14:01:15 -0700 Subject: [PATCH 11/14] Reflect SubMesh updates from MFEM v4.6. --- src/component_topology_handler.cpp | 10 ++++++---- src/topology_handler.cpp | 15 ++++++++++----- test/inputs/test_topol.2d.yml | 4 ++-- test/inputs/test_topol.3d.yml | 3 --- test/test_topol.cpp | 3 +++ 5 files changed, 21 insertions(+), 14 deletions(-) diff --git a/src/component_topology_handler.cpp b/src/component_topology_handler.cpp index 1af0ebfd..660af315 100644 --- a/src/component_topology_handler.cpp +++ b/src/component_topology_handler.cpp @@ -179,9 +179,9 @@ void ComponentTopologyHandler::SetupReferencePorts() for (int p = 0; p < port_list.size(); p++) { // Read hdf5 files. - std::string filename = config.GetRequiredOptionFromDict("file", port_list[p]); + std::string filename = config.GetOptionFromDict("file", "", port_list[p]); - if (FileExists(filename)) + if ((filename != "") && FileExists(filename)) ReadPortDatasFromFile(filename); else BuildPortDataFromInput(port_list[p]); @@ -320,9 +320,11 @@ void ComponentTopologyHandler::ReadPortsFromFile(const std::string filename) int attr_offset = 0; for (int c = 0; c < num_comp; c++) attr_offset = max(attr_offset, components[c]->bdr_attributes.Max()); - // Also does not overlap with global boundary attributes. assert(bdr_attributes.Size() > 0); - attr_offset = max(attr_offset, bdr_attributes.Max()); + // Also does not overlap with global boundary attributes. + // In order to keep consistencity with SubMesh, + // We leave bdr_attributes.Max()+1 empty. + attr_offset = max(attr_offset, bdr_attributes.Max()+1); attr_offset += 1; for (int p = 0; p < num_ports; p++) diff --git a/src/topology_handler.cpp b/src/topology_handler.cpp index 60954a0a..2c3ba34d 100644 --- a/src/topology_handler.cpp +++ b/src/topology_handler.cpp @@ -306,6 +306,9 @@ void SubMeshTopologyHandler::BuildSubMeshBoundary2D(const Mesh& pm, SubMesh& sm, MFEM_ASSERT(pm.Dimension() == 2, "Support only 2-dimension meshes!"); MFEM_ASSERT(sm.Dimension() == 2, "Support only 2-dimension meshes!"); + // NOTE(kevin): MFEM v4.6 SubMesh uses this for generated boundary attributes. + const int generated_battr = pmesh->bdr_attributes.Max() + 1; + // Array parent_face_map = submesh.GetParentFaceIDMap(); if (parent_face_map == NULL) parent_face_map = new Array(BuildFaceMap2D(pm, sm)); @@ -332,7 +335,7 @@ void SubMeshTopologyHandler::BuildSubMeshBoundary2D(const Mesh& pm, SubMesh& sm, // This case happens when a domain is extracted, but the root parent // mesh didn't have a boundary element on the surface that defined // it's boundary. It still creates a valid mesh, so we allow it. - sm.GetBdrElement(k)->SetAttribute(SubMesh::GENERATED_ATTRIBUTE); + sm.GetBdrElement(k)->SetAttribute(generated_battr); } } @@ -350,13 +353,15 @@ void SubMeshTopologyHandler::BuildInterfaceInfos() port_infos.SetSize(0); // interface attribute starts after the parent mesh boundary attributes. - int if_attr = pmesh->bdr_attributes.Max() + 1; + // NOTE(kevin): MFEM v4.6 SubMesh uses this for generated boundary attributes. + const int generated_battr = pmesh->bdr_attributes.Max() + 1; + int if_attr = pmesh->bdr_attributes.Max() + 2; for (int i = 0; i < numSub; i++) { for (int ib = 0; ib < meshes[i]->GetNBE(); ib++) { - if (meshes[i]->GetBdrAttribute(ib) != SubMesh::GENERATED_ATTRIBUTE) continue; + if (meshes[i]->GetBdrAttribute(ib) != generated_battr) continue; int parent_face_i = (*parent_face_map[i])[meshes[i]->GetBdrFace(ib)]; // Loop over each subdomain, each boundary element, to find the match. @@ -367,12 +372,12 @@ void SubMeshTopologyHandler::BuildInterfaceInfos() int parent_face_j = (*parent_face_map[j])[meshes[j]->GetBdrFace(jb)]; if (parent_face_i != parent_face_j) continue; - assert(meshes[j]->GetBdrAttribute(jb) == SubMesh::GENERATED_ATTRIBUTE); + assert(meshes[j]->GetBdrAttribute(jb) == generated_battr); if (interface_attributes[i][j] <= 0) { interface_attributes[i][j] = if_attr; interface_map[i][j] = port_infos.Size(); - // NOTE: for SubMehs, component boundary attribute is simply SubMesh::GENERATED_ATTRIBUTE, + // NOTE: for SubMesh, component boundary attribute is simply generated_battr, // which is not informative at all. port_infos.Append(PortInfo({.PortAttr = if_attr, .Mesh1 = i, .Mesh2 = j, diff --git a/test/inputs/test_topol.2d.yml b/test/inputs/test_topol.2d.yml index e1ca3bdf..cc15436f 100644 --- a/test/inputs/test_topol.2d.yml +++ b/test/inputs/test_topol.2d.yml @@ -9,7 +9,7 @@ mesh: file: "meshes/test.2x2.mesh" ports: - name: "port1" - file: "port1.h5" + file: "port1.2d.h5" comp1: name: "unit" attr: 3 @@ -18,7 +18,7 @@ mesh: attr: 1 comp2_configuration: [0.0, 2.0, 0.0, 0.0, 0.0, 0.0] - name: "port2" - file: "port2.h5" + file: "port2.2d.h5" comp1: name: "unit" attr: 2 diff --git a/test/inputs/test_topol.3d.yml b/test/inputs/test_topol.3d.yml index 2bf8bd2c..50d7082b 100644 --- a/test/inputs/test_topol.3d.yml +++ b/test/inputs/test_topol.3d.yml @@ -9,7 +9,6 @@ mesh: file: "meshes/test.1x1x1.hex.mesh" ports: - name: "port1" - file: "port1.h5" comp1: name: "cube" attr: 3 @@ -18,7 +17,6 @@ mesh: attr: 5 comp2_configuration: [0.5, 0.0, 0.0, 0.0, 0.0, 0.0] - name: "port2" - file: "port2.h5" comp1: name: "cube" attr: 4 @@ -27,7 +25,6 @@ mesh: attr: 2 comp2_configuration: [0.0, 0.5, 0.0, 0.0, 0.0, 0.0] - name: "port3" - file: "port3.h5" comp1: name: "cube" attr: 6 diff --git a/test/test_topol.cpp b/test/test_topol.cpp index c87dfe35..ca0924ce 100644 --- a/test/test_topol.cpp +++ b/test/test_topol.cpp @@ -62,6 +62,9 @@ TEST(PortReadWrite_test, Test_topol) { config = InputParser("inputs/test_topol.3d.yml"); config.dict_["mesh"]["component-wise"]["write_ports"] = true; + config.dict_["mesh"]["component-wise"]["ports"][0]["file"] = "port1.3d.h5"; + config.dict_["mesh"]["component-wise"]["ports"][1]["file"] = "port2.3d.h5"; + config.dict_["mesh"]["component-wise"]["ports"][2]["file"] = "port3.3d.h5"; printf("Generate\n"); ComponentTopologyHandler *old = new ComponentTopologyHandler(); From 0e7cf6709ef4be036f0ac1c83d836d7ca95f91d0 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Mon, 18 Mar 2024 12:18:35 -0700 Subject: [PATCH 12/14] visualization includes flow field. --- include/advdiff_solver.hpp | 9 ++++++ src/advdiff_solver.cpp | 56 ++++++++++++++++++++++++++++++++++++-- 2 files changed, 63 insertions(+), 2 deletions(-) diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index 303e7662..deea7f43 100644 --- a/include/advdiff_solver.hpp +++ b/include/advdiff_solver.hpp @@ -31,6 +31,13 @@ friend class ParameterizedProblem; bool save_flow = false; std::string flow_file = ""; + /* grid functions for visualizaing flow field */ + /* NOTE(kevin): this will be set up at SaveVisualization. */ + Array flow_fes; + Array flow_visual; + FiniteElementSpace *global_flow_fes = NULL; + Array global_flow_visual; + public: AdvDiffSolver(); @@ -47,6 +54,8 @@ friend class ParameterizedProblem; void SetParameterizedProblem(ParameterizedProblem *problem) override; + void SaveVisualization() override; + protected: void SetMUMPSSolver() override; diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 61bd190d..57ca4245 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -11,7 +11,7 @@ using namespace std; using namespace mfem; AdvDiffSolver::AdvDiffSolver() - : PoissonSolver() + : PoissonSolver(), flow_visual(0), flow_fes(0), global_flow_visual(0) { // ConvectionIntegrator does not support L2 space. assert(!full_dg); @@ -30,6 +30,9 @@ AdvDiffSolver::AdvDiffSolver() AdvDiffSolver::~AdvDiffSolver() { DeletePointers(flow_coeffs); + if (!stokes_solver) DeletePointers(flow_visual); + DeletePointers(global_flow_visual); + delete global_flow_fes; delete stokes_solver; } @@ -173,6 +176,56 @@ void AdvDiffSolver::SetParameterizedProblem(ParameterizedProblem *problem) PoissonSolver::SetParameterizedProblem(problem); } +void AdvDiffSolver::SaveVisualization() +{ + if (!save_visual) return; + + assert(paraviewColls.Size() > 0); + for (int k = 0; k < paraviewColls.Size(); k++) + assert(paraviewColls[k]); + + flow_visual.SetSize(numSub); + flow_visual = NULL; + + const int stokes_numvar = (stokes_solver) ? stokes_solver->GetNumVar() : 0; + if (!stokes_solver) + { + flow_fes.SetSize(numSub); + for (int m = 0; m < numSub; m++) + { + flow_fes[m] = new FiniteElementSpace(meshes[m], fec[0], dim); + flow_visual[m] = new GridFunction(flow_fes[m]); + } + } + + for (int m = 0; m < numSub; m++) + { + if (stokes_solver) + flow_visual[m] = stokes_solver->GetGridFunction(m * stokes_numvar); + else + { + assert(flow_coeffs[m]); + flow_visual[m]->ProjectCoefficient(*flow_coeffs[m]); + } + } + + if (unified_paraview) + { + assert(pmesh); + global_flow_fes = new FiniteElementSpace(pmesh, fec[0], dim); + global_flow_visual.SetSize(1); + global_flow_visual[0] = new GridFunction(global_flow_fes); + topol_handler->TransferToGlobal(flow_visual, global_flow_visual, 1); + + paraviewColls[0]->RegisterField("flow", global_flow_visual[0]); + } + else + for (int m = 0; m < numSub; m++) + paraviewColls[m]->RegisterField("flow", flow_visual[m]); + + MultiBlockSolver::SaveVisualization(); +} + void AdvDiffSolver::SetMUMPSSolver() { assert(globalMat_hypre); @@ -188,7 +241,6 @@ void AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem) stokes_solver = new StokesSolver; stokes_solver->InitVariables(); - stokes_solver->InitVisualization(); if (load_flow && FileExists(flow_file)) stokes_solver->LoadSolution(flow_file); From 272c85c2accb51876df5b17fd8cc205d466dd009 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Mon, 18 Mar 2024 15:28:43 -0700 Subject: [PATCH 13/14] save stokes flow solution --- include/multiblock_solver.hpp | 2 ++ src/advdiff_solver.cpp | 1 + src/multiblock_solver.cpp | 8 +++++++- 3 files changed, 10 insertions(+), 1 deletion(-) diff --git a/include/multiblock_solver.hpp b/include/multiblock_solver.hpp index b1f6d2e0..8c9b98b2 100644 --- a/include/multiblock_solver.hpp +++ b/include/multiblock_solver.hpp @@ -137,6 +137,8 @@ friend class ParameterizedProblem; BlockVector* GetSolution() { return U; } BlockVector* GetSolutionCopy() { return new BlockVector(*U); } + void SetSolutionSaveMode(const bool save_sol_); + void GetVariableVector(const int &var_idx, BlockVector &global, BlockVector &var); void SetVariableVector(const int &var_idx, BlockVector &var, BlockVector &global); diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 57ca4245..242d87f5 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -241,6 +241,7 @@ void AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem) stokes_solver = new StokesSolver; stokes_solver->InitVariables(); + stokes_solver->SetSolutionSaveMode(save_flow); if (load_flow && FileExists(flow_file)) stokes_solver->LoadSolution(flow_file); diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index daa06802..c0f07553 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -125,7 +125,13 @@ void MultiBlockSolver::ParseInputs() train_mode = SetTrainMode(); // save solution if single run. - save_sol = config.GetOption("save_solution/enabled", false); + SetSolutionSaveMode(config.GetOption("save_solution/enabled", false)); +} + +void MultiBlockSolver::SetSolutionSaveMode(const bool save_sol_) +{ + // save solution if single run. + save_sol = save_sol_; if (save_sol) { // Default file path if no input file name is provided. From 73c32ea0007a9815b863fc8ce787afb547ce8bd9 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Mon, 18 Mar 2024 15:43:30 -0700 Subject: [PATCH 14/14] save flow field only when it is newly computed. --- src/advdiff_solver.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 242d87f5..c8524d93 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -243,8 +243,12 @@ void AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem) stokes_solver->InitVariables(); stokes_solver->SetSolutionSaveMode(save_flow); + bool flow_loaded = false; if (load_flow && FileExists(flow_file)) + { stokes_solver->LoadSolution(flow_file); + flow_loaded = true; + } else { stokes_solver->SetParameterizedProblem(flow_problem); @@ -255,7 +259,7 @@ void AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem) stokes_solver->Solve(); } - if (save_flow) + if (save_flow && (!flow_loaded)) stokes_solver->SaveSolution(flow_file); DeletePointers(flow_coeffs);