From 6eba5ee198664ccac4d31ec94cb8d136f7925d37 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 31 Jan 2024 14:38:31 -0800 Subject: [PATCH 01/31] Add initial eigenvalue example --- CMakeLists.txt | 2 + examples/prom/laplace_eigenproblem.cpp | 526 +++++++++++++++++++++++++ 2 files changed, 528 insertions(+) create mode 100644 examples/prom/laplace_eigenproblem.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 8356f3e92..0be1b5c97 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -141,6 +141,7 @@ if (ENABLE_EXAMPLES) de_parametric_maxwell_greedy maxwell_local_rom_greedy grad_div_global_rom + laplace_eigenproblem dg_advection nonlinear_elasticity heat_conduction @@ -166,6 +167,7 @@ if (ENABLE_EXAMPLES) prom prom prom + prom dmd dmd dmd diff --git a/examples/prom/laplace_eigenproblem.cpp b/examples/prom/laplace_eigenproblem.cpp new file mode 100644 index 000000000..3abb23343 --- /dev/null +++ b/examples/prom/laplace_eigenproblem.cpp @@ -0,0 +1,526 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2024, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// libROM MFEM Example: Laplace Eigenproblem (adapted from ex11p.cpp) +// +// Compile with: make heat_conduction +// +// ================================================================================= +// +// Description: This example solves a time dependent nonlinear heat equation +// problem of the form du/dt = C(u), with a non-linear diffusion +// operator C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u. +// +// The example demonstrates the use of nonlinear operators (the +// class ConductionOperator defining C(u)), as well as their +// implicit time integration. Note that implementing the method +// ConductionOperator::ImplicitSolve is the only requirement for +// high-order implicit (SDIRK) time integration. Optional saving +// with ADIOS2 (adios2.readthedocs.io) is also illustrated. + +#include "mfem.hpp" +#include "linalg/BasisGenerator.h" +#include "linalg/BasisReader.h" +#include "linalg/Vector.h" +#include "mfem/Utilities.hpp" +#include +#include +#include + +using namespace std; +using namespace mfem; + +double v_initial(const Vector &x); +double kappa; + +int main(int argc, char *argv[]) +{ + // 1. Initialize MPI. + int num_procs, myid; + MPI_Init(&argc, &argv); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + + // 2. Parse command-line options. + const char *mesh_file = "../data/star.mesh"; + int ser_ref_levels = 2; + int par_ref_levels = 1; + int order = 1; + int nev = 5; + int seed = 75; + bool slu_solver = false; + bool sp_solver = false; + bool cpardiso_solver = false; + + double alpha = 1.0; + bool visualization = true; + bool visit = false; + int vis_steps = 5; + bool fom = false; + bool offline = false; + bool merge = false; + bool online = false; + int id = 0; + int nsets = 0; + + int precision = 8; + cout.precision(precision); + + OptionsParser args(argc, argv); + args.AddOption(&mesh_file, "-m", "--mesh", + "Mesh file to use."); + args.AddOption(&ser_ref_levels, "-rs", "--refine-serial", + "Number of times to refine the mesh uniformly in serial."); + args.AddOption(&par_ref_levels, "-rp", "--refine-parallel", + "Number of times to refine the mesh uniformly in parallel."); + args.AddOption(&order, "-o", "--order", + "Order (degree) of the finite elements."); + args.AddOption(&nev, "-n", "--num-eigs", + "Number of desired eigenmodes."); + args.AddOption(&seed, "-s", "--seed", + "Random seed used to initialize LOBPCG."); + args.AddOption(&alpha, "-a", "--alpha", + "Alpha coefficient."); + args.AddOption(&id, "-id", "--id", "Parametric id"); + args.AddOption(&nsets, "-ns", "--nset", "Number of parametric snapshot sets"); + args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", + "--no-visualization", + "Enable or disable GLVis visualization."); + args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit", + "--no-visit-datafiles", + "Save data files for VisIt (visit.llnl.gov) visualization."); + args.AddOption(&vis_steps, "-vs", "--visualization-steps", + "Visualize every n-th timestep."); + args.AddOption(&fom, "-fom", "--fom", "-no-fom", "--no-fom", + "Enable or disable the fom phase."); + args.AddOption(&offline, "-offline", "--offline", "-no-offline", "--no-offline", + "Enable or disable the offline phase."); + args.AddOption(&online, "-online", "--online", "-no-online", "--no-online", + "Enable or disable the online phase."); + args.AddOption(&merge, "-merge", "--merge", "-no-merge", "--no-merge", + "Enable or disable the merge phase."); + args.Parse(); + if (!args.Good()) + { + args.PrintUsage(cout); + MPI_Finalize(); + return 1; + } + + if (myid == 0) + { + args.PrintOptions(cout); + } + + kappa = alpha * M_PI; + + if (fom) + { + MFEM_VERIFY(fom && !offline && !online + && !merge, "everything must be turned off if fom is used."); + } + else + { + bool check = (offline && !merge && !online) || (!offline && merge && !online) + || (!offline && !merge && online); + MFEM_VERIFY(check, "only one of offline, merge, or online must be true!"); + } + + // 3. Read the serial mesh from the given mesh file on all processors. We can + // handle triangular, quadrilateral, tetrahedral and hexahedral meshes + // with the same code. + Mesh *mesh = new Mesh(mesh_file, 1, 1); + int dim = mesh->Dimension(); + + // 4. Refine the mesh in serial to increase the resolution. In this example + // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is + // a command-line parameter. + for (int lev = 0; lev < ser_ref_levels; lev++) + { + mesh->UniformRefinement(); + } + + // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine + // this mesh further in parallel to increase the resolution. Once the + // parallel mesh is defined, the serial mesh can be deleted. + ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh); + delete mesh; + for (int lev = 0; lev < par_ref_levels; lev++) + { + pmesh->UniformRefinement(); + } + + // 6. Define a parallel finite element space on the parallel mesh. Here we + // use continuous Lagrange finite elements of the specified order. If + // order < 1, we instead use an isoparametric/isogeometric space. + FiniteElementCollection *fec; + if (order > 0) + { + fec = new H1_FECollection(order, dim); + } + else if (pmesh->GetNodes()) + { + fec = pmesh->GetNodes()->OwnFEC(); + } + else + { + fec = new H1_FECollection(order = 1, dim); + } + ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec); + HYPRE_BigInt size = fespace->GlobalTrueVSize(); + if (myid == 0) + { + cout << "Number of unknowns: " << size << endl; + } + + // 7. Initiate ROM related variables + int max_num_snapshots = 100; + bool update_right_SV = false; + bool isIncremental = false; + const std::string basisName = "basis"; + const std::string basisFileName = basisName + std::to_string(id); + CAROM::Options* options; + CAROM::BasisGenerator *generator; + StopWatch solveTimer, assembleTimer, mergeTimer; + + // 8. Set BasisGenerator if offline + if (offline) + { + options = new CAROM::Options(fespace->GetTrueVSize(), nev, nev, + update_right_SV); + generator = new CAROM::BasisGenerator(*options, isIncremental, basisFileName); + } + + // 9. The merge phase + if (merge) + { + mergeTimer.Start(); + options = new CAROM::Options(fespace->GetTrueVSize(), max_num_snapshots, nev, + update_right_SV); + generator = new CAROM::BasisGenerator(*options, isIncremental, basisName); + for (int paramID=0; paramIDloadSamples(snapshot_filename, "snapshot"); + } + generator->endSamples(); // save the merged basis file + mergeTimer.Stop(); + if (myid == 0) + { + printf("Elapsed time for merging and building ROM basis: %e second\n", + mergeTimer.RealTime()); + } + delete generator; + delete options; + MPI_Finalize(); + return 0; + } + + // 10. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite + // element space. The first corresponds to the Laplacian operator -Delta, + // while the second is a simple mass matrix needed on the right hand side + // of the generalized eigenvalue problem below. The boundary conditions + // are implemented by elimination with special values on the diagonal to + // shift the Dirichlet eigenvalues out of the computational range. After + // serial and parallel assembly we extract the corresponding parallel + // matrices A and M. + ConstantCoefficient one(1.0); + + ParGridFunction u_gf(fespace); + FunctionCoefficient u_0(v_initial); + u_gf.ProjectCoefficient(u_0); + + Vector u; + u_gf.GetTrueDofs(u); + + Array ess_bdr; + if (pmesh->bdr_attributes.Size()) + { + ess_bdr.SetSize(pmesh->bdr_attributes.Max()); + ess_bdr = 1; + } + + ParBilinearForm *a = new ParBilinearForm(fespace); + a->AddDomainIntegrator(new DiffusionIntegrator(u_0)); + if (pmesh->bdr_attributes.Size() == 0) + { + // Add a mass term if the mesh has no boundary, e.g. periodic mesh or + // closed surface. + a->AddDomainIntegrator(new MassIntegrator(one)); + } + a->Assemble(); + a->EliminateEssentialBCDiag(ess_bdr, 1.0); + a->Finalize(); + + ParBilinearForm *m = new ParBilinearForm(fespace); + m->AddDomainIntegrator(new MassIntegrator(one)); + m->Assemble(); + // shift the eigenvalue corresponding to eliminated dofs to a large value + m->EliminateEssentialBCDiag(ess_bdr, numeric_limits::min()); + m->Finalize(); + + HypreParMatrix *A = a->ParallelAssemble(); + HypreParMatrix *M = m->ParallelAssemble(); + + delete a; + delete m; + + ParGridFunction x(fespace); + HypreLOBPCG *lobpcg; + Array eigenvalues; + DenseMatrix evect; + + // 11. The offline phase + if(fom || offline) + { + // 12. Define and configure the LOBPCG eigensolver and the BoomerAMG + // preconditioner for A to be used within the solver. Set the matrices + // which define the generalized eigenproblem A x = lambda M x. + Solver * precond = NULL; + if (!slu_solver && !sp_solver && !cpardiso_solver) + { + HypreBoomerAMG * amg = new HypreBoomerAMG(*A); + amg->SetPrintLevel(1); + precond = amg; + } + else + { + // TODO: preconditioners using MFEM_SUPERLU, STRUMPACK, CPARDISO + } + lobpcg = new HypreLOBPCG(MPI_COMM_WORLD); + lobpcg->SetNumModes(nev); + lobpcg->SetRandomSeed(seed); + lobpcg->SetPreconditioner(*precond); + lobpcg->SetMaxIter(200); + lobpcg->SetTol(1e-8); + lobpcg->SetPrecondUsageMode(1); + lobpcg->SetPrintLevel(1); + lobpcg->SetMassMatrix(*M); + lobpcg->SetOperator(*A); + + // 13. Compute the eigenmodes and extract the array of eigenvalues. Define a + // parallel grid function to represent each of the eigenmodes returned by + // the solver. + solveTimer.Start(); + lobpcg->Solve(); + solveTimer.Stop(); + lobpcg->GetEigenvalues(eigenvalues); + + // 14. take and write snapshots for ROM + if (offline) + { + for (int i=0; iGetEigenvector(i); + generator->takeSample(x.GetData(), (double)i, 0.01); + } + generator->writeSnapshot(); + + delete generator; + delete options; + } + + delete precond; + } + + // 15. The online phase + if (online) { + // 16. read the reduced basis + assembleTimer.Start(); + CAROM::BasisReader reader(basisName); + + Vector ev; + const CAROM::Matrix* spatialbasis = reader.getSpatialBasis(0.0); + const int numRowRB = spatialbasis->numRows(); + const int numColumnRB = spatialbasis->numColumns(); + if (myid == 0) printf("spatial basis dimension is %d x %d\n", numRowRB, + numColumnRB); + + // 17. form ROM operator + CAROM::Matrix invReducedA; + ComputeCtAB(*A, *spatialbasis, *spatialbasis, invReducedA); + DenseMatrix *A_mat = new DenseMatrix(invReducedA.numRows(), + invReducedA.numColumns()); + A_mat->Set(1, invReducedA.getData()); + + CAROM::Matrix invReducedM; + ComputeCtAB(*M, *spatialbasis, *spatialbasis, invReducedM); + DenseMatrix *M_mat = new DenseMatrix(invReducedM.numRows(), + invReducedM.numColumns()); + M_mat->Set(1, invReducedM.getData()); + + assembleTimer.Stop(); + + // 18. solve ROM + solveTimer.Start(); + // (Q^T A Q) c = \lamba (Q^T M Q) c + A_mat->Eigenvalues(*M_mat, ev, evect); + solveTimer.Stop(); + + if (myid == 0) + { + eigenvalues = Array(ev.GetData(), ev.Size()); + for (int i = 0; i < ev.Size(); i++) + { + std::cout << "Eigenvalue " << i << ": = " << eigenvalues[i] << "\n"; + } + } + + delete spatialbasis; + + delete A_mat; + delete M_mat; + } + + // 19. Save the refined mesh and the modes in parallel. This output can be + // viewed later using GLVis: "glvis -np -m mesh -g mode". + { + ostringstream mesh_name, mode_name; + mesh_name << "laplace_eigenproblem-mesh." << setfill('0') << setw(6) << myid; + + ofstream mesh_ofs(mesh_name.str().c_str()); + mesh_ofs.precision(8); + pmesh->Print(mesh_ofs); + + for (int i=0; iGetEigenvector(i); + } else { + // for online, eigenvectors are stored in evect matrix + Vector ev; + evect.GetRow(i, ev); + x = ev; + } + + mode_name << "mode_" << setfill('0') << setw(2) << i << "." + << setfill('0') << setw(6) << myid; + + ofstream mode_ofs(mode_name.str().c_str()); + mode_ofs.precision(8); + x.Save(mode_ofs); + mode_name.str(""); + } + } + + VisItDataCollection visit_dc("LaplaceEigenproblem", pmesh); + if (visit) + { + visit_dc.RegisterField("v", &u_gf); + for (int i=0; iGetEigenvector(i); + } else { + // for online, eigenvectors are stored in evect matrix + Vector ev; + evect.GetRow(i, ev); + x = ev; + } + visit_dc.RegisterField("x" + std::to_string(i), &x); + } + visit_dc.SetCycle(0); + visit_dc.SetTime(0.0); + visit_dc.Save(); + } + + socketstream sout; + if (visualization) + { + char vishost[] = "localhost"; + int visport = 19916; + sout.open(vishost, visport); + sout << "parallel " << num_procs << " " << myid << endl; + int good = sout.good(), all_good; + MPI_Allreduce(&good, &all_good, 1, MPI_INT, MPI_MIN, pmesh->GetComm()); + if (!all_good) + { + sout.close(); + visualization = false; + if (myid == 0) + { + cout << "Unable to connect to GLVis server at " + << vishost << ':' << visport << endl; + cout << "GLVis visualization disabled.\n"; + } + } + else + { + for (int i=0; i " << flush; + cin >> c; + } + MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD); + + if (c != 'c') + { + break; + } + } + sout.close(); + } + } + + // 20. Free the used memory. + if (fom || offline) + { + delete lobpcg; + } + delete M; + delete A; + + delete fespace; + if (order > 0) + { + delete fec; + } + delete pmesh; + + MPI_Finalize(); + + return 0; +} + +double v_initial(const Vector &x) +{ + return 1.0 + cos(kappa*x(1))*sin(kappa*x(0)); +} From cc46d068cc81955b2cde77e721224bccdc47d1d1 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Fri, 2 Feb 2024 09:04:33 -0800 Subject: [PATCH 02/31] Rename example and update description --- CMakeLists.txt | 2 +- ...pp => laplace_eigenproblem_global_rom.cpp} | 33 ++++++++++++------- 2 files changed, 23 insertions(+), 12 deletions(-) rename examples/prom/{laplace_eigenproblem.cpp => laplace_eigenproblem_global_rom.cpp} (93%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0be1b5c97..ca7d8499b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -141,7 +141,7 @@ if (ENABLE_EXAMPLES) de_parametric_maxwell_greedy maxwell_local_rom_greedy grad_div_global_rom - laplace_eigenproblem + laplace_eigenproblem_global_rom dg_advection nonlinear_elasticity heat_conduction diff --git a/examples/prom/laplace_eigenproblem.cpp b/examples/prom/laplace_eigenproblem_global_rom.cpp similarity index 93% rename from examples/prom/laplace_eigenproblem.cpp rename to examples/prom/laplace_eigenproblem_global_rom.cpp index 3abb23343..e14727e04 100644 --- a/examples/prom/laplace_eigenproblem.cpp +++ b/examples/prom/laplace_eigenproblem_global_rom.cpp @@ -10,20 +10,31 @@ // libROM MFEM Example: Laplace Eigenproblem (adapted from ex11p.cpp) // -// Compile with: make heat_conduction -// // ================================================================================= // -// Description: This example solves a time dependent nonlinear heat equation -// problem of the form du/dt = C(u), with a non-linear diffusion -// operator C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u. +// Description: This example code demonstrates the use of MFEM and libROM to +// define a simple projection-based reduced order model of the +// eigenvalue problem -Delta ((1 + alpha v) u) = lambda u with homogeneous +// Dirichlet boundary conditions where alpha is a scalar ROM parameter +// controlling the frequency of v. +// +// We compute a number of the lowest eigenmodes by discretizing +// the Laplacian and Mass operators using a FE space of the +// specified order, or an isoparametric/isogeometric space if +// order < 1 (quadratic for quadratic curvilinear mesh, NURBS for +// NURBS mesh, etc.) +// +// Offline phase: laplace_eigenproblem_global_rom -offline -id 0 -a 0.8 +// laplace_eigenproblem_global_rom -offline -id 1 -a 0.9 +// laplace_eigenproblem_global_rom -offline -id 2 -a 1.1 +// laplace_eigenproblem_global_rom -offline -id 3 -a 1.2 +// +// Merge phase: laplace_eigenproblem_global_rom -merge -ns 4 +// +// FOM run (for error calculation): +// laplace_eigenproblem_global_rom -fom -f 1.0 // -// The example demonstrates the use of nonlinear operators (the -// class ConductionOperator defining C(u)), as well as their -// implicit time integration. Note that implementing the method -// ConductionOperator::ImplicitSolve is the only requirement for -// high-order implicit (SDIRK) time integration. Optional saving -// with ADIOS2 (adios2.readthedocs.io) is also illustrated. +// Online phase: laplace_eigenproblem_global_rom -online -f 1.0 #include "mfem.hpp" #include "linalg/BasisGenerator.h" From 9f9eff3c1324e50df34f7029269e2feafe8921e0 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Fri, 2 Feb 2024 09:39:17 -0800 Subject: [PATCH 03/31] Add error calculation with FOM --- .../prom/laplace_eigenproblem_global_rom.cpp | 63 ++++++++++++++++--- 1 file changed, 54 insertions(+), 9 deletions(-) diff --git a/examples/prom/laplace_eigenproblem_global_rom.cpp b/examples/prom/laplace_eigenproblem_global_rom.cpp index e14727e04..ce3e79d1a 100644 --- a/examples/prom/laplace_eigenproblem_global_rom.cpp +++ b/examples/prom/laplace_eigenproblem_global_rom.cpp @@ -313,7 +313,7 @@ int main(int argc, char *argv[]) lobpcg->SetMaxIter(200); lobpcg->SetTol(1e-8); lobpcg->SetPrecondUsageMode(1); - lobpcg->SetPrintLevel(1); + lobpcg->SetPrintLevel(0); lobpcg->SetMassMatrix(*M); lobpcg->SetOperator(*A); @@ -326,19 +326,22 @@ int main(int argc, char *argv[]) lobpcg->GetEigenvalues(eigenvalues); // 14. take and write snapshots for ROM - if (offline) + for (int i = 0; i < nev; i++) { - for (int i=0; iGetEigenvector(i); generator->takeSample(x.GetData(), (double)i, 0.01); } - generator->writeSnapshot(); + } + if (offline) + { + generator->writeSnapshot(); delete generator; delete options; } @@ -383,7 +386,7 @@ int main(int argc, char *argv[]) if (myid == 0) { eigenvalues = Array(ev.GetData(), ev.Size()); - for (int i = 0; i < ev.Size(); i++) + for (int i = 0; i < nev; i++) { std::cout << "Eigenvalue " << i << ": = " << eigenvalues[i] << "\n"; } @@ -395,6 +398,41 @@ int main(int argc, char *argv[]) delete M_mat; } + ostringstream sol_ev_name, sol_ev_name_fom; + if (fom || offline) + { + sol_ev_name << "sol_eigenvalues_fom." << setfill('0') << setw(6) << myid; + } + if (online) + { + sol_ev_name << "sol_eigenvalues." << setfill('0') << setw(6) << myid; + sol_ev_name_fom << "sol_eigenvalues_fom." << setfill('0') << setw(6) << myid; + } + + if (online) + { + // Initialize FOM solution + Vector ev_fom(nev); + + ifstream fom_file; + fom_file.open(sol_ev_name_fom.str().c_str()); + ev_fom.Load(fom_file, ev_fom.Size()); + fom_file.close(); + + Vector diff_ev(nev); + for (int i = 0; i < nev; i++) + { + diff_ev[i] = ev_fom[i] - eigenvalues[i]; + double ev_diff_norm = sqrt(diff_ev[i] * diff_ev[i]); + double ev_fom_norm = sqrt(ev_fom[i] * ev_fom[i]); + if (myid == 0) + { + std::cout << "Relative error of ROM solution for eigenvalue " << i << " = " << + ev_diff_norm / ev_fom_norm << std::endl; + } + } + } + // 19. Save the refined mesh and the modes in parallel. This output can be // viewed later using GLVis: "glvis -np -m mesh -g mode". { @@ -425,6 +463,13 @@ int main(int argc, char *argv[]) x.Save(mode_ofs); mode_name.str(""); } + + ofstream sol_ev_ofs(sol_ev_name.str().c_str()); + sol_ev_ofs.precision(16); + for (int i = 0; i < nev; ++i) + { + sol_ev_ofs << eigenvalues[i] << std::endl; + } } VisItDataCollection visit_dc("LaplaceEigenproblem", pmesh); From 03ec16f2d62a74233f555f89f7f4a807d3edf627 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Fri, 2 Feb 2024 10:06:27 -0800 Subject: [PATCH 04/31] Add options for energy fraction and reduced dim --- .../prom/laplace_eigenproblem_global_rom.cpp | 29 ++++++++++++++----- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/examples/prom/laplace_eigenproblem_global_rom.cpp b/examples/prom/laplace_eigenproblem_global_rom.cpp index ce3e79d1a..cc45d92b3 100644 --- a/examples/prom/laplace_eigenproblem_global_rom.cpp +++ b/examples/prom/laplace_eigenproblem_global_rom.cpp @@ -80,6 +80,8 @@ int main(int argc, char *argv[]) bool online = false; int id = 0; int nsets = 0; + double ef = 0.9999; + int rdim = -1; int precision = 8; cout.precision(precision); @@ -117,6 +119,10 @@ int main(int argc, char *argv[]) "Enable or disable the online phase."); args.AddOption(&merge, "-merge", "--merge", "-no-merge", "--no-merge", "Enable or disable the merge phase."); + args.AddOption(&ef, "-ef", "--energy_fraction", + "Energy fraction."); + args.AddOption(&rdim, "-rdim", "--rdim", + "Reduced dimension."); args.Parse(); if (!args.Good()) { @@ -356,7 +362,16 @@ int main(int argc, char *argv[]) CAROM::BasisReader reader(basisName); Vector ev; - const CAROM::Matrix* spatialbasis = reader.getSpatialBasis(0.0); + const CAROM::Matrix *spatialbasis; + if (rdim != -1) + { + spatialbasis = reader.getSpatialBasis(0.0, rdim); + } + else + { + spatialbasis = reader.getSpatialBasis(0.0, ef); + } + const int numRowRB = spatialbasis->numRows(); const int numColumnRB = spatialbasis->numColumns(); if (myid == 0) printf("spatial basis dimension is %d x %d\n", numRowRB, @@ -386,7 +401,7 @@ int main(int argc, char *argv[]) if (myid == 0) { eigenvalues = Array(ev.GetData(), ev.Size()); - for (int i = 0; i < nev; i++) + for (int i = 0; i < ev.Size() && i < nev; i++) { std::cout << "Eigenvalue " << i << ": = " << eigenvalues[i] << "\n"; } @@ -420,7 +435,7 @@ int main(int argc, char *argv[]) fom_file.close(); Vector diff_ev(nev); - for (int i = 0; i < nev; i++) + for (int i = 0; i < eigenvalues.Size() && i < nev; i++) { diff_ev[i] = ev_fom[i] - eigenvalues[i]; double ev_diff_norm = sqrt(diff_ev[i] * diff_ev[i]); @@ -443,7 +458,7 @@ int main(int argc, char *argv[]) mesh_ofs.precision(8); pmesh->Print(mesh_ofs); - for (int i=0; i Date: Fri, 2 Feb 2024 15:47:55 -0800 Subject: [PATCH 05/31] Update eigenproblem problem setup --- .../prom/laplace_eigenproblem_global_rom.cpp | 72 ++++++++++++++----- 1 file changed, 55 insertions(+), 17 deletions(-) diff --git a/examples/prom/laplace_eigenproblem_global_rom.cpp b/examples/prom/laplace_eigenproblem_global_rom.cpp index cc45d92b3..68109d45d 100644 --- a/examples/prom/laplace_eigenproblem_global_rom.cpp +++ b/examples/prom/laplace_eigenproblem_global_rom.cpp @@ -48,8 +48,11 @@ using namespace std; using namespace mfem; -double v_initial(const Vector &x); +double Conductivity(const Vector &x); +double Potential(const Vector &x); +int problem = 1; double kappa; +Vector bb_min, bb_max; int main(int argc, char *argv[]) { @@ -89,6 +92,8 @@ int main(int argc, char *argv[]) OptionsParser args(argc, argv); args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use."); + args.AddOption(&problem, "-p", "--problem", + "Problem setup to use. See options in Conductivity and Potential functions."); args.AddOption(&ser_ref_levels, "-rs", "--refine-serial", "Number of times to refine the mesh uniformly in serial."); args.AddOption(&par_ref_levels, "-rp", "--refine-parallel", @@ -251,12 +256,8 @@ int main(int argc, char *argv[]) // matrices A and M. ConstantCoefficient one(1.0); - ParGridFunction u_gf(fespace); - FunctionCoefficient u_0(v_initial); - u_gf.ProjectCoefficient(u_0); - - Vector u; - u_gf.GetTrueDofs(u); + FunctionCoefficient kappa_0(Conductivity); + FunctionCoefficient v_0(Potential); Array ess_bdr; if (pmesh->bdr_attributes.Size()) @@ -266,13 +267,8 @@ int main(int argc, char *argv[]) } ParBilinearForm *a = new ParBilinearForm(fespace); - a->AddDomainIntegrator(new DiffusionIntegrator(u_0)); - if (pmesh->bdr_attributes.Size() == 0) - { - // Add a mass term if the mesh has no boundary, e.g. periodic mesh or - // closed surface. - a->AddDomainIntegrator(new MassIntegrator(one)); - } + a->AddDomainIntegrator(new DiffusionIntegrator(kappa_0)); + a->AddDomainIntegrator(new MassIntegrator(v_0)); a->Assemble(); a->EliminateEssentialBCDiag(ess_bdr, 1.0); a->Finalize(); @@ -490,7 +486,6 @@ int main(int argc, char *argv[]) VisItDataCollection visit_dc("LaplaceEigenproblem", pmesh); if (visit) { - visit_dc.RegisterField("v", &u_gf); for (int i = 0; i < nev && i < eigenvalues.Size(); i++) { if (fom || offline) { @@ -591,7 +586,50 @@ int main(int argc, char *argv[]) return 0; } -double v_initial(const Vector &x) +double Conductivity(const Vector &x) +{ + int dim = x.Size(); + + switch (problem) + { + case 1: + return 1.0; + case 2: + return 1.0 + cos(kappa * x(1)) * sin(kappa * x(0)); + case 3: + case 4: + case 5: + return 1.0; + } + return 0.0; +} + +double Potential(const Vector &x) { - return 1.0 + cos(kappa*x(1))*sin(kappa*x(0)); + int dim = x.Size(); + + // map x to the reference [-1,1] domain + Vector X(dim); + for (int i = 0; i < dim; i++) + { + double center = (bb_min[i] + bb_max[i]) * 0.5; + X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]); + } + + Vector center(dim); + center = kappa; + + switch (problem) + { + case 1: + case 2: + return 0.0; + case 3: + return 1.0; + case 4: + return std::exp(X.DistanceSquaredTo(center) / 0.01); + case 5: + return -std::exp(X.DistanceSquaredTo(center) / 0.01); + } + return 0.0; } From 92a702c98b94b76fe7fda3f9efd2cb11b5206e0f Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Tue, 6 Feb 2024 13:23:38 -0800 Subject: [PATCH 06/31] Fix center and kappa for problem 4 and 5 --- examples/prom/laplace_eigenproblem_global_rom.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/prom/laplace_eigenproblem_global_rom.cpp b/examples/prom/laplace_eigenproblem_global_rom.cpp index 68109d45d..261537f4a 100644 --- a/examples/prom/laplace_eigenproblem_global_rom.cpp +++ b/examples/prom/laplace_eigenproblem_global_rom.cpp @@ -168,6 +168,7 @@ int main(int argc, char *argv[]) { mesh->UniformRefinement(); } + mesh->GetBoundingBox(bb_min, bb_max, max(order, 1)); // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine // this mesh further in parallel to increase the resolution. Once the @@ -612,12 +613,12 @@ double Potential(const Vector &x) Vector X(dim); for (int i = 0; i < dim; i++) { - double center = (bb_min[i] + bb_max[i]) * 0.5; - X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]); + double mesh_center = (bb_min[i] + bb_max[i]) * 0.5; + X(i) = 2.0 * (x(i) - mesh_center) / (bb_max[i] - bb_min[i]); } - Vector center(dim); - center = kappa; + Vector center(dim); // center of gaussian for problem 4 and 5 + center = kappa / M_PI; switch (problem) { From 66f86a3744585f61203bc21e21bc9c964924d09f Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 7 Feb 2024 09:08:50 -0800 Subject: [PATCH 07/31] Fix saving eigenvectors to visit --- examples/prom/laplace_eigenproblem_global_rom.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/examples/prom/laplace_eigenproblem_global_rom.cpp b/examples/prom/laplace_eigenproblem_global_rom.cpp index 261537f4a..1a5e40b4b 100644 --- a/examples/prom/laplace_eigenproblem_global_rom.cpp +++ b/examples/prom/laplace_eigenproblem_global_rom.cpp @@ -487,6 +487,7 @@ int main(int argc, char *argv[]) VisItDataCollection visit_dc("LaplaceEigenproblem", pmesh); if (visit) { + std::vector visit_evs; for (int i = 0; i < nev && i < eigenvalues.Size(); i++) { if (fom || offline) { @@ -498,11 +499,16 @@ int main(int argc, char *argv[]) evect.GetRow(i, ev); x = ev; } - visit_dc.RegisterField("x" + std::to_string(i), &x); + visit_evs.push_back(new ParGridFunction(x)); + visit_dc.RegisterField("x" + std::to_string(i), visit_evs.back()); } visit_dc.SetCycle(0); visit_dc.SetTime(0.0); visit_dc.Save(); + for (size_t i = 0; i < visit_evs.size(); i++) + { + delete visit_evs[i]; + } } socketstream sout; From e0285fec8e3efd48bbdc18c6db7c1798b73e39a1 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 7 Feb 2024 09:21:16 -0800 Subject: [PATCH 08/31] Rename example to elliptic eigenvalue --- CMakeLists.txt | 2 +- ...p => elliptic_eigenproblem_global_rom.cpp} | 20 +++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) rename examples/prom/{laplace_eigenproblem_global_rom.cpp => elliptic_eigenproblem_global_rom.cpp} (96%) diff --git a/CMakeLists.txt b/CMakeLists.txt index ca7d8499b..6b39d1955 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -141,7 +141,7 @@ if (ENABLE_EXAMPLES) de_parametric_maxwell_greedy maxwell_local_rom_greedy grad_div_global_rom - laplace_eigenproblem_global_rom + elliptic_eigenproblem_global_rom dg_advection nonlinear_elasticity heat_conduction diff --git a/examples/prom/laplace_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp similarity index 96% rename from examples/prom/laplace_eigenproblem_global_rom.cpp rename to examples/prom/elliptic_eigenproblem_global_rom.cpp index 1a5e40b4b..dd1d6c011 100644 --- a/examples/prom/laplace_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -8,7 +8,7 @@ * *****************************************************************************/ -// libROM MFEM Example: Laplace Eigenproblem (adapted from ex11p.cpp) +// libROM MFEM Example: Elliptic Eigenproblem (adapted from ex11p.cpp) // // ================================================================================= // @@ -24,17 +24,17 @@ // order < 1 (quadratic for quadratic curvilinear mesh, NURBS for // NURBS mesh, etc.) // -// Offline phase: laplace_eigenproblem_global_rom -offline -id 0 -a 0.8 -// laplace_eigenproblem_global_rom -offline -id 1 -a 0.9 -// laplace_eigenproblem_global_rom -offline -id 2 -a 1.1 -// laplace_eigenproblem_global_rom -offline -id 3 -a 1.2 +// Offline phase: elliptic_eigenproblem_global_rom -offline -p 2 -id 0 -a 0.8 +// elliptic_eigenproblem_global_rom -offline -p 2 -id 1 -a 0.9 +// elliptic_eigenproblem_global_rom -offline -p 2 -id 2 -a 1.1 +// elliptic_eigenproblem_global_rom -offline -p 2 -id 3 -a 1.2 // -// Merge phase: laplace_eigenproblem_global_rom -merge -ns 4 +// Merge phase: elliptic_eigenproblem_global_rom -merge -p 2 -ns 4 // // FOM run (for error calculation): -// laplace_eigenproblem_global_rom -fom -f 1.0 +// elliptic_eigenproblem_global_rom -fom -p 2 -f 1.0 // -// Online phase: laplace_eigenproblem_global_rom -online -f 1.0 +// Online phase: elliptic_eigenproblem_global_rom -online -p 2 -f 1.0 #include "mfem.hpp" #include "linalg/BasisGenerator.h" @@ -449,7 +449,7 @@ int main(int argc, char *argv[]) // viewed later using GLVis: "glvis -np -m mesh -g mode". { ostringstream mesh_name, mode_name; - mesh_name << "laplace_eigenproblem-mesh." << setfill('0') << setw(6) << myid; + mesh_name << "elliptic_eigenproblem-mesh." << setfill('0') << setw(6) << myid; ofstream mesh_ofs(mesh_name.str().c_str()); mesh_ofs.precision(8); @@ -484,7 +484,7 @@ int main(int argc, char *argv[]) } } - VisItDataCollection visit_dc("LaplaceEigenproblem", pmesh); + VisItDataCollection visit_dc("EllipticEigenproblem", pmesh); if (visit) { std::vector visit_evs; From ed00f76317aa632216dd4d089ee96f5920cde0f9 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 7 Feb 2024 10:53:52 -0800 Subject: [PATCH 09/31] Update case 4 and 5 --- examples/prom/elliptic_eigenproblem_global_rom.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index dd1d6c011..769385b97 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -634,9 +634,9 @@ double Potential(const Vector &x) case 3: return 1.0; case 4: - return std::exp(X.DistanceSquaredTo(center) / 0.01); + return std::exp(-X.DistanceSquaredTo(center) / 0.01); case 5: - return -std::exp(X.DistanceSquaredTo(center) / 0.01); + return -std::exp(-X.DistanceSquaredTo(center) / 0.01); } return 0.0; } From 108ce8cb75a557f8c301f3b5304862f0a159b11f Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Thu, 8 Feb 2024 08:10:29 -0800 Subject: [PATCH 10/31] Output initial condition and adjust center relative to mesh --- .../prom/elliptic_eigenproblem_global_rom.cpp | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 769385b97..63a9a5912 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -260,6 +260,13 @@ int main(int argc, char *argv[]) FunctionCoefficient kappa_0(Conductivity); FunctionCoefficient v_0(Potential); + // Project initial conductivity and potential to visualize initial condition + ParGridFunction c_gf(fespace); + c_gf.ProjectCoefficient(kappa_0); + + ParGridFunction p_gf(fespace); + p_gf.ProjectCoefficient(v_0); + Array ess_bdr; if (pmesh->bdr_attributes.Size()) { @@ -487,6 +494,8 @@ int main(int argc, char *argv[]) VisItDataCollection visit_dc("EllipticEigenproblem", pmesh); if (visit) { + visit_dc.RegisterField("InitialConductivity", &c_gf); + visit_dc.RegisterField("InitialPotential", &p_gf); std::vector visit_evs; for (int i = 0; i < nev && i < eigenvalues.Size(); i++) { @@ -617,15 +626,15 @@ double Potential(const Vector &x) // map x to the reference [-1,1] domain Vector X(dim); + Vector center(dim); // center of gaussian for problem 4 and 5 + center = kappa / M_PI; for (int i = 0; i < dim; i++) { double mesh_center = (bb_min[i] + bb_max[i]) * 0.5; X(i) = 2.0 * (x(i) - mesh_center) / (bb_max[i] - bb_min[i]); + center(i) = 2.0 * (center(i) - mesh_center) / (bb_max[i] - bb_min[i]); } - Vector center(dim); // center of gaussian for problem 4 and 5 - center = kappa / M_PI; - switch (problem) { case 1: From ca421b2cda1fd04906f2cbdb5419e33614766275 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Thu, 8 Feb 2024 12:59:42 -0800 Subject: [PATCH 11/31] Add case 6 and use square mesh by default --- .../prom/elliptic_eigenproblem_global_rom.cpp | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 63a9a5912..db825ae67 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -63,7 +63,7 @@ int main(int argc, char *argv[]) MPI_Comm_rank(MPI_COMM_WORLD, &myid); // 2. Parse command-line options. - const char *mesh_file = "../data/star.mesh"; + const char *mesh_file = ""; int ser_ref_levels = 2; int par_ref_levels = 1; int order = 1; @@ -158,7 +158,15 @@ int main(int argc, char *argv[]) // 3. Read the serial mesh from the given mesh file on all processors. We can // handle triangular, quadrilateral, tetrahedral and hexahedral meshes // with the same code. - Mesh *mesh = new Mesh(mesh_file, 1, 1); + Mesh *mesh; + if (mesh_file == "") + { + mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL)); + } + else + { + mesh = new Mesh(mesh_file, 1, 1); + } int dim = mesh->Dimension(); // 4. Refine the mesh in serial to increase the resolution. In this example @@ -615,6 +623,7 @@ double Conductivity(const Vector &x) case 3: case 4: case 5: + case 6: return 1.0; } return 0.0; @@ -643,9 +652,14 @@ double Potential(const Vector &x) case 3: return 1.0; case 4: - return std::exp(-X.DistanceSquaredTo(center) / 0.01); case 5: return -std::exp(-X.DistanceSquaredTo(center) / 0.01); + case 6: + const double D = 1.0; + Vector neg_center(center); + neg_center.Neg(); + return -D * (std::exp(-X.DistanceSquaredTo(center) / 0.01) + std::exp( + -X.DistanceSquaredTo(neg_center) / 0.01)); } return 0.0; } From 6328f6bd06f00def75254b1fcb073a95681268b0 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Tue, 13 Feb 2024 14:01:20 -0800 Subject: [PATCH 12/31] Update mesh bounds and mapping --- .../prom/elliptic_eigenproblem_global_rom.cpp | 34 +++++++++++++------ 1 file changed, 23 insertions(+), 11 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index db825ae67..9edaf501c 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -161,7 +161,8 @@ int main(int argc, char *argv[]) Mesh *mesh; if (mesh_file == "") { - mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL)); + mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false, + 100.0, 100.0)); } else { @@ -633,17 +634,30 @@ double Potential(const Vector &x) { int dim = x.Size(); - // map x to the reference [-1,1] domain Vector X(dim); Vector center(dim); // center of gaussian for problem 4 and 5 center = kappa / M_PI; + + Vector neg_center(center); + neg_center.Neg(); + + // map x to the reference [-1,1] domain for (int i = 0; i < dim; i++) { - double mesh_center = (bb_min[i] + bb_max[i]) * 0.5; - X(i) = 2.0 * (x(i) - mesh_center) / (bb_max[i] - bb_min[i]); - center(i) = 2.0 * (center(i) - mesh_center) / (bb_max[i] - bb_min[i]); + X(i) = bb_min[i] + (bb_max[i] - bb_min[i]) * ((x(i) - bb_min[i]) / + (bb_max[i] - bb_min[i])); + + // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) + center(i) = bb_min[i] + (center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * 0.5); + neg_center(i) = bb_min[i] + (neg_center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * + 0.5); } + // amplitude of gaussians for problems 4-6 + const double D = 1.0; + // width of gaussians for problems 4-6 + const double c = 10.0; + switch (problem) { case 1: @@ -652,14 +666,12 @@ double Potential(const Vector &x) case 3: return 1.0; case 4: + return D * std::exp(-X.DistanceSquaredTo(center) / c); case 5: - return -std::exp(-X.DistanceSquaredTo(center) / 0.01); + return -D * std::exp(-X.DistanceSquaredTo(center) / c); case 6: - const double D = 1.0; - Vector neg_center(center); - neg_center.Neg(); - return -D * (std::exp(-X.DistanceSquaredTo(center) / 0.01) + std::exp( - -X.DistanceSquaredTo(neg_center) / 0.01)); + return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp( + -X.DistanceSquaredTo(neg_center) / c)); } return 0.0; } From 95ed0264c6da79772dc84188102469cfc31f74e9 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 14 Feb 2024 07:55:37 -0800 Subject: [PATCH 13/31] Add additional case to vary center around point --- .../prom/elliptic_eigenproblem_global_rom.cpp | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 9edaf501c..16c6097a4 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -625,6 +625,7 @@ double Conductivity(const Vector &x) case 4: case 5: case 6: + case 7: return 1.0; } return 0.0; @@ -672,6 +673,26 @@ double Potential(const Vector &x) case 6: return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp( -X.DistanceSquaredTo(neg_center) / c)); + case 7: + center = 0.25; + center(0) += 0.05 * cos(2.0*kappa); + center(1) += 0.05 * sin(2.0*kappa); + + neg_center = center; + neg_center.Neg(); + + center.Print(); + neg_center.Print(); + + for (int i = 0; i < dim; i++) + { + // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) + center(i) = bb_min[i] + (center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * 0.5); + neg_center(i) = bb_min[i] + (neg_center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * + 0.5); + } + return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp( + -X.DistanceSquaredTo(neg_center) / c)); } return 0.0; } From f76f825045353990851de73adbf377dff1f12bcb Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 14 Feb 2024 11:58:07 -0800 Subject: [PATCH 14/31] Transform mesh to be centered around origin --- .../prom/elliptic_eigenproblem_global_rom.cpp | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 16c6097a4..ed7665f0d 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -161,8 +161,22 @@ int main(int argc, char *argv[]) Mesh *mesh; if (mesh_file == "") { + double x_max = 20.0; + double y_max = 20.0; mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false, - 100.0, 100.0)); + x_max, y_max)); + mesh->GetBoundingBox(bb_min, bb_max, max(order, 1)); + + // shift mesh around origin (from [0, bb_max] -> [-bb_max/2, bb_max/2]) + auto *mesh_transform = +[](const Vector &v_orig, Vector &v_new) -> void { + v_new = v_orig; + // shift mesh vertices by (bb_max[i] / 2) + v_new(0) -= 0.5*bb_max[0]; + v_new(1) -= 0.5*bb_max[1]; + }; + + // performs the transform - bb_min/bb_max are updated again below + mesh->Transform(mesh_transform); } else { @@ -681,9 +695,6 @@ double Potential(const Vector &x) neg_center = center; neg_center.Neg(); - center.Print(); - neg_center.Print(); - for (int i = 0; i < dim; i++) { // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) From ba1e591869c189cc0b6347a839eb3fbc9cce2917 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 14 Feb 2024 12:05:15 -0800 Subject: [PATCH 15/31] Add modified cases for limiting domain --- .../prom/elliptic_eigenproblem_global_rom.cpp | 62 ++++++++++++++++++- 1 file changed, 60 insertions(+), 2 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index ed7665f0d..7c829f468 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -640,6 +640,8 @@ double Conductivity(const Vector &x) case 5: case 6: case 7: + case 8: + case 9: return 1.0; } return 0.0; @@ -669,9 +671,28 @@ double Potential(const Vector &x) } // amplitude of gaussians for problems 4-6 - const double D = 1.0; + const double D = 100.0; // width of gaussians for problems 4-6 - const double c = 10.0; + const double c = 0.05; + + const double min_d = 2.0 * sqrt(2.0); + + auto check_domain = [&](const Vector &t, const double limit) { + // helper function to check if t is within limit percentage relative to the center of the mesh + for (int i = 0; i < dim; i++) + { + double domain_limit = limit * (bb_max[i] - bb_min[i]); + double mesh_center = 0.5 * (bb_max[i] + bb_min[i]); + + // check that t is within the limit relative to the center of the mesh + if (t(i) - mesh_center > 0.5 * domain_limit) + { + std::cerr << "Error: value of t exceeds domain limit: t = " << t( + i) << ", limit = " << 0.5 * domain_limit << "\n"; + exit(-1); + } + } + }; switch (problem) { @@ -704,6 +725,43 @@ double Potential(const Vector &x) } return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp( -X.DistanceSquaredTo(neg_center) / c)); + case 8: + // Similar to case 6, but t is restricted to inner 20% of the domain + // in this case, the radius of the gaussian is (0.05*min_d)^2 where + // min_d is the lower bound of the atom distance over time + + // verify t is within inner 20% of domain (centered around mesh origin) + check_domain(center, 0.2); + + return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d, + 2)) + std::exp( + -X.DistanceSquaredTo(neg_center) / std::pow(c * min_d, 2))); + case 9: + // Similar to case 7, but t is restricted to inner 20% of the domain and t is defined as: + // t = (1.5 + 0.5*cos(2*k), 1.5 + 0.5*sin(2*k)) where k = alpha * PI, alpha is the input parameter given with the -a option. + // The radius of the gaussian follows case 8: (0.05*min_d)^2 + + center = 0.15; + center(0) += 0.05 * cos(2.0*kappa); + center(1) += 0.05 * sin(2.0*kappa); + + neg_center = center; + neg_center.Neg(); + + for (int i = 0; i < dim; i++) + { + // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) + center(i) = bb_min[i] + (center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * 0.5); + neg_center(i) = bb_min[i] + (neg_center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * + 0.5); + } + + // verify t is within inner 20% of domain (centered around mesh origin) + check_domain(center, 0.2); + + return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d, + 2)) + std::exp( + -X.DistanceSquaredTo(neg_center) / std::pow(c * min_d, 2))); } return 0.0; } From dff666c3032d0add6b624eb3ca087736987efec0 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 14 Feb 2024 12:21:59 -0800 Subject: [PATCH 16/31] Refactor parameter mapping used in potential setup --- .../prom/elliptic_eigenproblem_global_rom.cpp | 44 +++++++++++-------- 1 file changed, 26 insertions(+), 18 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 7c829f468..4d2c56891 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -658,18 +658,6 @@ double Potential(const Vector &x) Vector neg_center(center); neg_center.Neg(); - // map x to the reference [-1,1] domain - for (int i = 0; i < dim; i++) - { - X(i) = bb_min[i] + (bb_max[i] - bb_min[i]) * ((x(i) - bb_min[i]) / - (bb_max[i] - bb_min[i])); - - // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) - center(i) = bb_min[i] + (center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * 0.5); - neg_center(i) = bb_min[i] + (neg_center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * - 0.5); - } - // amplitude of gaussians for problems 4-6 const double D = 100.0; // width of gaussians for problems 4-6 @@ -694,6 +682,28 @@ double Potential(const Vector &x) } }; + auto map_fraction_to_mesh = [](const double &bb_min, const double &bb_max, + const double &fraction) -> double { + // helper function to map a fractional value from [-1, 1] to [bb_min, bb_max] + CAROM_VERIFY(fraction <= 1.0 && fraction >= -1.0); + return bb_min + (fraction + 1.0) * ((bb_max - bb_min) * 0.5); + }; + + auto map_mesh_to_fraction = [](const double &bb_min, const double &bb_max, + const double &value) -> double { + // helper function to map a value from the mesh range [bb_min, bb_max] to [-1, 1] + CAROM_VERIFY(value <= bb_max && value >= bb_min); + return -1.0 + (value - bb_min) * ((2.0) / (bb_max - bb_min)); + }; + + X = x; + for (int i = 0; i < dim; i++) + { + // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) + center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], center(i)); + neg_center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], neg_center(i)); + } + switch (problem) { case 1: @@ -719,9 +729,8 @@ double Potential(const Vector &x) for (int i = 0; i < dim; i++) { // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) - center(i) = bb_min[i] + (center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * 0.5); - neg_center(i) = bb_min[i] + (neg_center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * - 0.5); + center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], center(i)); + neg_center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], neg_center(i)); } return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp( -X.DistanceSquaredTo(neg_center) / c)); @@ -751,9 +760,8 @@ double Potential(const Vector &x) for (int i = 0; i < dim; i++) { // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) - center(i) = bb_min[i] + (center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * 0.5); - neg_center(i) = bb_min[i] + (neg_center(i) + 1.0) * ((bb_max[i] - bb_min[i]) * - 0.5); + center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], center(i)); + neg_center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], neg_center(i)); } // verify t is within inner 20% of domain (centered around mesh origin) From c089e15848d555bc2da0d29c47fdb2bc849a1051 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 14 Feb 2024 13:33:53 -0800 Subject: [PATCH 17/31] Update case 9 setup to use absolute mesh location --- examples/prom/elliptic_eigenproblem_global_rom.cpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 4d2c56891..b0b447fd9 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -750,9 +750,14 @@ double Potential(const Vector &x) // t = (1.5 + 0.5*cos(2*k), 1.5 + 0.5*sin(2*k)) where k = alpha * PI, alpha is the input parameter given with the -a option. // The radius of the gaussian follows case 8: (0.05*min_d)^2 - center = 0.15; - center(0) += 0.05 * cos(2.0*kappa); - center(1) += 0.05 * sin(2.0*kappa); + // Sets the center to vary around +/- 1.5 (absolute location on the mesh) + center = 1.5; + center(0) += 0.5 * cos(2.0*kappa); + center(1) += 0.5 * sin(2.0*kappa); + + // map the absolute location back to a fraction of the mesh domain + center(0) = map_mesh_to_fraction(bb_min[0], bb_max[0], center(0)); + center(1) = map_mesh_to_fraction(bb_min[1], bb_max[1], center(1)); neg_center = center; neg_center.Neg(); From 5d51359bf41d881ce3639a6653b23cc5cafb82a5 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 14 Feb 2024 17:53:21 -0800 Subject: [PATCH 18/31] Add eigenvector error calculation --- .../prom/elliptic_eigenproblem_global_rom.cpp | 114 +++++++++++++----- 1 file changed, 87 insertions(+), 27 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index b0b447fd9..3992162fe 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -434,6 +434,18 @@ int main(int argc, char *argv[]) } } + DenseMatrix tmp(evect); + evect = DenseMatrix(nev, numRowRB); + for (int i = 0; i < ev.Size() && i < nev; i++) + { + Vector evector; + tmp.GetRow(i, evector); + CAROM::Vector evector_carom(evector.GetData(), evector.Size(), false, false); + CAROM::Vector *ev_carom = spatialbasis->mult(evector_carom); + evect.SetRow(i, ev_carom->getData()); + delete ev_carom; + } + delete spatialbasis; delete A_mat; @@ -451,30 +463,6 @@ int main(int argc, char *argv[]) sol_ev_name_fom << "sol_eigenvalues_fom." << setfill('0') << setw(6) << myid; } - if (online) - { - // Initialize FOM solution - Vector ev_fom(nev); - - ifstream fom_file; - fom_file.open(sol_ev_name_fom.str().c_str()); - ev_fom.Load(fom_file, ev_fom.Size()); - fom_file.close(); - - Vector diff_ev(nev); - for (int i = 0; i < eigenvalues.Size() && i < nev; i++) - { - diff_ev[i] = ev_fom[i] - eigenvalues[i]; - double ev_diff_norm = sqrt(diff_ev[i] * diff_ev[i]); - double ev_fom_norm = sqrt(ev_fom[i] * ev_fom[i]); - if (myid == 0) - { - std::cout << "Relative error of ROM solution for eigenvalue " << i << " = " << - ev_diff_norm / ev_fom_norm << std::endl; - } - } - } - // 19. Save the refined mesh and the modes in parallel. This output can be // viewed later using GLVis: "glvis -np -m mesh -g mode". { @@ -485,6 +473,15 @@ int main(int argc, char *argv[]) mesh_ofs.precision(8); pmesh->Print(mesh_ofs); + std::string mode_prefix = "mode_"; + if (online) + { + mode_prefix += "rom_"; + } else if (fom) + { + mode_prefix += "fom_"; + } + for (int i=0; i < nev && i < eigenvalues.Size(); i++) { if (fom || offline) { @@ -497,12 +494,19 @@ int main(int argc, char *argv[]) x = ev; } - mode_name << "mode_" << setfill('0') << setw(2) << i << "." + mode_name << mode_prefix << setfill('0') << setw(2) << i << "." << setfill('0') << setw(6) << myid; ofstream mode_ofs(mode_name.str().c_str()); - mode_ofs.precision(8); - x.Save(mode_ofs); + mode_ofs.precision(16); + + // TODO: issue using .Load() function if file written with .Save()? + //x.Save(mode_ofs); + for (int j = 0; j < x.Size(); j++) + { + mode_ofs << x[j] << "\n"; + } + mode_ofs.close(); mode_name.str(""); } @@ -514,6 +518,62 @@ int main(int argc, char *argv[]) } } + if (online) + { + // Initialize FOM solution + Vector ev_fom(nev); + + ifstream fom_file; + fom_file.open(sol_ev_name_fom.str().c_str()); + ev_fom.Load(fom_file, ev_fom.Size()); + fom_file.close(); + + Vector diff_ev(nev); + for (int i = 0; i < eigenvalues.Size() && i < nev; i++) + { + diff_ev[i] = ev_fom[i] - eigenvalues[i]; + double ev_diff_norm = sqrt(diff_ev[i] * diff_ev[i]); + double ev_fom_norm = sqrt(ev_fom[i] * ev_fom[i]); + if (myid == 0) + { + std::cout << "Relative error of ROM solution for eigenvalue " << i << " = " << + ev_diff_norm / ev_fom_norm << std::endl; + } + } + + // Calculate errors of eigenvectors + ostringstream mode_name, mode_name_fom; + Vector mode_rom(evect.NumCols()); + Vector mode_fom(evect.NumCols()); + for (int i = 0; i < eigenvalues.Size() && i < nev; i++) + { + mode_name << "mode_rom_" << setfill('0') << setw(2) << i << "." + << setfill('0') << setw(6) << myid; + mode_name_fom << "mode_fom_" << setfill('0') << setw(2) << i << "." + << setfill('0') << setw(6) << myid; + + ifstream mode_rom_ifs(mode_name.str().c_str()); + mode_rom_ifs.precision(16); + mode_rom.Load(mode_rom_ifs, evect.NumCols()); + mode_rom_ifs.close(); + + ifstream mode_fom_ifs(mode_name_fom.str().c_str()); + mode_fom_ifs.precision(16); + mode_fom.Load(mode_fom_ifs, evect.NumCols()); + mode_fom_ifs.close(); + + const double fomNorm = sqrt(InnerProduct(mode_fom, mode_fom)); + mode_fom -= mode_rom; + const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom)); + if (myid == 0) std::cout << "Relative l2 error of ROM eigenvector " << i << + " = " << diffNorm / + fomNorm << std::endl; + + mode_name.str(""); + mode_name_fom.str(""); + } + } + VisItDataCollection visit_dc("EllipticEigenproblem", pmesh); if (visit) { From 8bb701132272ac0fce6165cfcb65f5f417f68887 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Thu, 15 Feb 2024 08:17:00 -0800 Subject: [PATCH 19/31] Add timing info printout --- .../prom/elliptic_eigenproblem_global_rom.cpp | 24 +++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 3992162fe..0b1a60a86 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -278,6 +278,7 @@ int main(int argc, char *argv[]) // shift the Dirichlet eigenvalues out of the computational range. After // serial and parallel assembly we extract the corresponding parallel // matrices A and M. + assembleTimer.Start(); ConstantCoefficient one(1.0); FunctionCoefficient kappa_0(Conductivity); @@ -314,6 +315,8 @@ int main(int argc, char *argv[]) HypreParMatrix *A = a->ParallelAssemble(); HypreParMatrix *M = m->ParallelAssemble(); + assembleTimer.Stop(); + delete a; delete m; @@ -665,7 +668,24 @@ int main(int argc, char *argv[]) } } - // 20. Free the used memory. + // 20. print timing info + if (myid == 0) + { + if(fom || offline) + { + printf("Elapsed time for assembling FOM: %e second\n", + assembleTimer.RealTime()); + printf("Elapsed time for solving FOM: %e second\n", solveTimer.RealTime()); + } + if(online) + { + printf("Elapsed time for assembling ROM: %e second\n", + assembleTimer.RealTime()); + printf("Elapsed time for solving ROM: %e second\n", solveTimer.RealTime()); + } + } + + // 21. Free the used memory. if (fom || offline) { delete lobpcg; @@ -733,7 +753,7 @@ double Potential(const Vector &x) double mesh_center = 0.5 * (bb_max[i] + bb_min[i]); // check that t is within the limit relative to the center of the mesh - if (t(i) - mesh_center > 0.5 * domain_limit) + if ((t(i) - mesh_center) - (0.5 * domain_limit) > 1.0e-14) { std::cerr << "Error: value of t exceeds domain limit: t = " << t( i) << ", limit = " << 0.5 * domain_limit << "\n"; From 11367db62c40ba1d602de406f5e973aca4e66bfe Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 28 Feb 2024 06:33:27 -0800 Subject: [PATCH 20/31] Add option to compile mfem with lapack enabled --- scripts/compile.sh | 12 +++++++++++- scripts/setup.sh | 4 ++-- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/scripts/compile.sh b/scripts/compile.sh index 6dcc2243f..299cfc14c 100755 --- a/scripts/compile.sh +++ b/scripts/compile.sh @@ -26,6 +26,7 @@ USE_MFEM="Off" UPDATE_LIBS=false INSTALL_SCALAPACK=false MFEM_USE_GSLIB="Off" +MFEM_USE_LAPACK="NO" cleanup_dependencies() { pushd . @@ -45,7 +46,7 @@ cleanup_dependencies() { # Get options -while getopts "ah:dh:gh:mh:t:uh:sh" o; +while getopts "ah:dh:gh:lh:mh:t:uh:sh" o; do case "${o}" in a) @@ -57,6 +58,9 @@ do g) MFEM_USE_GSLIB="On" ;; + l) + MFEM_USE_LAPACK="YES" + ;; m) USE_MFEM="On" ;; @@ -105,6 +109,12 @@ if [[ $MFEM_USE_GSLIB == "On" ]] && [[ ! -d "$HOME_DIR/dependencies/gslib" ]]; t fi export MFEM_USE_GSLIB +if [[ ${MFEM_USE_LAPACK} == "YES" ]] && [[ ${USE_MFEM} == "Off" ]]; then + echo "mfem lapack flag has no effect without mfem enabled" + exit 1 +fi +export MFEM_USE_LAPACK + SCRIPT_DIR=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd ) REPO_PREFIX=$( dirname $SCRIPT_DIR ) diff --git a/scripts/setup.sh b/scripts/setup.sh index 636052816..a3e108e86 100755 --- a/scripts/setup.sh +++ b/scripts/setup.sh @@ -115,7 +115,7 @@ if [[ $BUILD_TYPE == "Debug" ]]; then if [[ $UPDATE_LIBS == "true" ]]; then cd mfem_debug git pull - make pdebug -j 8 STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" + make pdebug -j 8 STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_LAPACK=${MFEM_USE_LAPACK:-"NO"} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" check_result $? mfem-debug-installation fi cd $LIB_DIR @@ -132,7 +132,7 @@ else # NOTE(kevin): v4.5.2-dev commit. This is the mfem version used by PyMFEM v4.5.2.0. # This version matching is required to support pylibROM-PyMFEM interface. git checkout 00b2a0705f647e17a1d4ffcb289adca503f28d42 - make parallel -j 8 STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" + make parallel -j 8 STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_LAPACK=${MFEM_USE_LAPACK:-"NO"} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" check_result $? mfem-parallel-installation fi cd $LIB_DIR From ae4589065d54c7c54653ab901f3b1a34b790e83e Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Tue, 5 Mar 2024 06:02:33 -0800 Subject: [PATCH 21/31] Update takeSample to use new function signature --- examples/prom/elliptic_eigenproblem_global_rom.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 0b1a60a86..d19cab61d 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -371,7 +371,7 @@ int main(int argc, char *argv[]) if (offline) { x = lobpcg->GetEigenvector(i); - generator->takeSample(x.GetData(), (double)i, 0.01); + generator->takeSample(x.GetData()); } } From 046ca9286f866518288c5c7177d439368999a780 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Tue, 5 Mar 2024 09:28:10 -0800 Subject: [PATCH 22/31] Add verbose option for solver --- examples/prom/elliptic_eigenproblem_global_rom.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index d19cab61d..9308b59eb 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -85,6 +85,7 @@ int main(int argc, char *argv[]) int nsets = 0; double ef = 0.9999; int rdim = -1; + int verbose_level = 0; int precision = 8; cout.precision(precision); @@ -128,6 +129,8 @@ int main(int argc, char *argv[]) "Energy fraction."); args.AddOption(&rdim, "-rdim", "--rdim", "Reduced dimension."); + args.AddOption(&verbose_level, "-v", "--verbose", + "Set the verbosity level of the LOBPCG solver and preconditioner. 0 is off."); args.Parse(); if (!args.Good()) { @@ -335,7 +338,7 @@ int main(int argc, char *argv[]) if (!slu_solver && !sp_solver && !cpardiso_solver) { HypreBoomerAMG * amg = new HypreBoomerAMG(*A); - amg->SetPrintLevel(1); + amg->SetPrintLevel(verbose_level); precond = amg; } else @@ -349,7 +352,7 @@ int main(int argc, char *argv[]) lobpcg->SetMaxIter(200); lobpcg->SetTol(1e-8); lobpcg->SetPrecondUsageMode(1); - lobpcg->SetPrintLevel(0); + lobpcg->SetPrintLevel(verbose_level); lobpcg->SetMassMatrix(*M); lobpcg->SetOperator(*A); From a84548569e2e5a081042f4dbedcf8b33e0d886ce Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Mon, 11 Mar 2024 10:54:35 -0700 Subject: [PATCH 23/31] Renaming mapping helper functions --- .../prom/elliptic_eigenproblem_global_rom.cpp | 27 ++++++++++--------- scripts/compile.sh | 6 ++--- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 9308b59eb..ff0b79647 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -164,6 +164,7 @@ int main(int argc, char *argv[]) Mesh *mesh; if (mesh_file == "") { + // square domain with side length 20 at origin (0,0) double x_max = 20.0; double y_max = 20.0; mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false, @@ -748,7 +749,7 @@ double Potential(const Vector &x) const double min_d = 2.0 * sqrt(2.0); - auto check_domain = [&](const Vector &t, const double limit) { + auto check_within_portion = [&](const Vector &t, const double limit) { // helper function to check if t is within limit percentage relative to the center of the mesh for (int i = 0; i < dim; i++) { @@ -765,14 +766,14 @@ double Potential(const Vector &x) } }; - auto map_fraction_to_mesh = [](const double &bb_min, const double &bb_max, + auto map_to_ref_mesh = [](const double &bb_min, const double &bb_max, const double &fraction) -> double { // helper function to map a fractional value from [-1, 1] to [bb_min, bb_max] CAROM_VERIFY(fraction <= 1.0 && fraction >= -1.0); return bb_min + (fraction + 1.0) * ((bb_max - bb_min) * 0.5); }; - auto map_mesh_to_fraction = [](const double &bb_min, const double &bb_max, + auto map_from_ref_mesh = [](const double &bb_min, const double &bb_max, const double &value) -> double { // helper function to map a value from the mesh range [bb_min, bb_max] to [-1, 1] CAROM_VERIFY(value <= bb_max && value >= bb_min); @@ -783,8 +784,8 @@ double Potential(const Vector &x) for (int i = 0; i < dim; i++) { // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) - center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], center(i)); - neg_center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], neg_center(i)); + center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], center(i)); + neg_center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i)); } switch (problem) @@ -812,8 +813,8 @@ double Potential(const Vector &x) for (int i = 0; i < dim; i++) { // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) - center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], center(i)); - neg_center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], neg_center(i)); + center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], center(i)); + neg_center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i)); } return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp( -X.DistanceSquaredTo(neg_center) / c)); @@ -823,7 +824,7 @@ double Potential(const Vector &x) // min_d is the lower bound of the atom distance over time // verify t is within inner 20% of domain (centered around mesh origin) - check_domain(center, 0.2); + check_within_portion(center, 0.2); return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d, 2)) + std::exp( @@ -839,8 +840,8 @@ double Potential(const Vector &x) center(1) += 0.5 * sin(2.0*kappa); // map the absolute location back to a fraction of the mesh domain - center(0) = map_mesh_to_fraction(bb_min[0], bb_max[0], center(0)); - center(1) = map_mesh_to_fraction(bb_min[1], bb_max[1], center(1)); + center(0) = map_from_ref_mesh(bb_min[0], bb_max[0], center(0)); + center(1) = map_from_ref_mesh(bb_min[1], bb_max[1], center(1)); neg_center = center; neg_center.Neg(); @@ -848,12 +849,12 @@ double Potential(const Vector &x) for (int i = 0; i < dim; i++) { // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) - center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], center(i)); - neg_center(i) = map_fraction_to_mesh(bb_min[i], bb_max[i], neg_center(i)); + center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], center(i)); + neg_center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i)); } // verify t is within inner 20% of domain (centered around mesh origin) - check_domain(center, 0.2); + check_within_portion(center, 0.2); return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d, 2)) + std::exp( diff --git a/scripts/compile.sh b/scripts/compile.sh index 299cfc14c..9e43a3a3f 100755 --- a/scripts/compile.sh +++ b/scripts/compile.sh @@ -26,7 +26,7 @@ USE_MFEM="Off" UPDATE_LIBS=false INSTALL_SCALAPACK=false MFEM_USE_GSLIB="Off" -MFEM_USE_LAPACK="NO" +MFEM_USE_LAPACK="Off" cleanup_dependencies() { pushd . @@ -59,7 +59,7 @@ do MFEM_USE_GSLIB="On" ;; l) - MFEM_USE_LAPACK="YES" + MFEM_USE_LAPACK="On" ;; m) USE_MFEM="On" @@ -109,7 +109,7 @@ if [[ $MFEM_USE_GSLIB == "On" ]] && [[ ! -d "$HOME_DIR/dependencies/gslib" ]]; t fi export MFEM_USE_GSLIB -if [[ ${MFEM_USE_LAPACK} == "YES" ]] && [[ ${USE_MFEM} == "Off" ]]; then +if [[ ${MFEM_USE_LAPACK} == "On" ]] && [[ ${USE_MFEM} == "Off" ]]; then echo "mfem lapack flag has no effect without mfem enabled" exit 1 fi From 32f9e4f609d706acad8815dfc5c4b588abd6f72e Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 13 Mar 2024 15:16:41 -0700 Subject: [PATCH 24/31] Add SuperLU, Strumpack, MKL CPardiso preconditioners --- .../prom/elliptic_eigenproblem_global_rom.cpp | 94 ++++++++++++++++++- 1 file changed, 89 insertions(+), 5 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index ff0b79647..33ebb13c8 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -131,12 +131,40 @@ int main(int argc, char *argv[]) "Reduced dimension."); args.AddOption(&verbose_level, "-v", "--verbose", "Set the verbosity level of the LOBPCG solver and preconditioner. 0 is off."); +#ifdef MFEM_USE_SUPERLU + args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu", + "--no-superlu", "Use the SuperLU Solver."); +#endif +#ifdef MFEM_USE_STRUMPACK + args.AddOption(&sp_solver, "-sp", "--strumpack", "-no-sp", + "--no-strumpack", "Use the STRUMPACK Solver."); +#endif +#ifdef MFEM_USE_MKL_CPARDISO + args.AddOption(&cpardiso_solver, "-cpardiso", "--cpardiso", "-no-cpardiso", + "--no-cpardiso", "Use the MKL CPardiso Solver."); +#endif args.Parse(); - if (!args.Good()) + if (slu_solver && sp_solver) { - args.PrintUsage(cout); - MPI_Finalize(); - return 1; + if (myid == 0) + std::cout << "WARNING: Both SuperLU and STRUMPACK have been selected," + << " please choose either one." << std::endl + << " Defaulting to SuperLU." << std::endl; + sp_solver = false; + } + // The command line options are also passed to the STRUMPACK + // solver. So do not exit if some options are not recognized. + if (!sp_solver) + { + if (!args.Good()) + { + if (myid == 0) + { + args.PrintUsage(cout); + } + MPI_Finalize(); + return 1; + } } if (myid == 0) @@ -319,6 +347,22 @@ int main(int argc, char *argv[]) HypreParMatrix *A = a->ParallelAssemble(); HypreParMatrix *M = m->ParallelAssemble(); +#if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK) + Operator * Arow = NULL; +#endif +#ifdef MFEM_USE_SUPERLU + if (slu_solver) + { + Arow = new SuperLURowLocMatrix(*A); + } +#endif +#ifdef MFEM_USE_STRUMPACK + if (sp_solver) + { + Arow = new STRUMPACKRowLocMatrix(*A); + } +#endif + assembleTimer.Stop(); delete a; @@ -344,8 +388,44 @@ int main(int argc, char *argv[]) } else { - // TODO: preconditioners using MFEM_SUPERLU, STRUMPACK, CPARDISO +#ifdef MFEM_USE_SUPERLU + if (slu_solver) + { + SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD); + superlu->SetPrintStatistics(verbose_level > 0 ? true : false); + superlu->SetSymmetricPattern(true); + superlu->SetColumnPermutation(superlu::PARMETIS); + superlu->SetOperator(*Arow); + precond = superlu; + } +#endif +#ifdef MFEM_USE_STRUMPACK + if (sp_solver) + { + STRUMPACKSolver * strumpack = new STRUMPACKSolver(MPI_COMM_WORLD, argc, argv); + strumpack->SetPrintFactorStatistics(true); + strumpack->SetPrintSolveStatistics(verbose_level > 0 ? true : false); + strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT); + strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS); + strumpack->SetMatching(strumpack::MatchingJob::NONE); + strumpack->SetCompression(strumpack::CompressionType::NONE); + strumpack->SetOperator(*Arow); + strumpack->SetFromCommandLine(); + precond = strumpack; + } +#endif +#ifdef MFEM_USE_MKL_CPARDISO + if (cpardiso_solver) + { + auto cpardiso = new CPardisoSolver(A->GetComm()); + cpardiso->SetMatrixType(CPardisoSolver::MatType::REAL_STRUCTURE_SYMMETRIC); + cpardiso->SetPrintLevel(verbose_level); + cpardiso->SetOperator(*A); + precond = cpardiso; + } +#endif } + lobpcg = new HypreLOBPCG(MPI_COMM_WORLD); lobpcg->SetNumModes(nev); lobpcg->SetRandomSeed(seed); @@ -704,6 +784,10 @@ int main(int argc, char *argv[]) } delete pmesh; +#if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK) + delete Arow; +#endif + MPI_Finalize(); return 0; From 90488de1ba34b8d5303f8a4f403f9e48a8ef4617 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Thu, 14 Mar 2024 07:05:38 -0700 Subject: [PATCH 25/31] Add SuperLU build to setup script --- scripts/setup.sh | 41 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/scripts/setup.sh b/scripts/setup.sh index a3e108e86..fb9374625 100755 --- a/scripts/setup.sh +++ b/scripts/setup.sh @@ -105,6 +105,43 @@ METIS_DIR=$LIB_DIR/parmetis-4.0.3 METIS_OPT=-I${METIS_DIR}/metis/include METIS_LIB="-L${METIS_DIR}/build/lib/libparmetis -lparmetis -L${METIS_DIR}/build/lib/libmetis -lmetis" +# Install distributed version of SuperLU +cd $LIB_DIR +if [ -z ${INSTALL_SUPERLU} ] && [ ! -d "superlu_dist" ]; then + wget https://github.com/xiaoyeli/superlu_dist/archive/refs/tags/v6.3.1.tar.gz + tar -zxvf v6.3.1.tar.gz + mv superlu_dist-6.3.1 superlu_dist + cd superlu_dist + mkdir build + cd build + export PARMETIS_ROOT=${METIS_DIR} + export PARMETIS_BUILD_DIR=${PARMETIS_ROOT}/build/lib/ + cmake -DCMAKE_INSTALL_PREFIX=./ \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_C_COMPILER=mpicc \ + -DCMAKE_CXX_COMPILER=mpicxx \ + -DCMAKE_Fortran_COMPILER=mpifort \ + -Denable_complex16=OFF \ + -Denable_examples=OFF \ + -DTPL_ENABLE_PARMETISLIB=ON \ + -DTPL_PARMETIS_INCLUDE_DIRS="${PARMETIS_ROOT}/include;${PARMETIS_ROOT}/metis/include" \ + -DTPL_PARMETIS_LIBRARIES="${PARMETIS_BUILD_DIR}/libparmetis/libparmetis.a;${PARMETIS_BUILD_DIR}/libmetis/libmetis.a" \ + -DTPL_ENABLE_COMBBLASLIB=OFF \ + -DUSE_XSDK_DEFAULTS="False" \ + -DBUILD_SHARED_LIBS=ON \ + -DBUILD_STATIC_LIBS=OFF \ + .. + check_result $? superlu-cmake + make -j 8 + check_result $? superlu-build + make install + check_result $? superlu-installation +fi + +SUPERLU_DIR=$LIB_DIR/superlu_dist +SUPERLU_OPT=-I${SUPERLU_DIR}/build/include +SUPERLU_LIB="-Wl,-rpath,${SUPERLU_DIR}/build/lib/ -L${SUPERLU_DIR}/build/lib/ -lsuperlu_dist -lblas" + # Install MFEM cd $LIB_DIR if [[ $BUILD_TYPE == "Debug" ]]; then @@ -115,7 +152,7 @@ if [[ $BUILD_TYPE == "Debug" ]]; then if [[ $UPDATE_LIBS == "true" ]]; then cd mfem_debug git pull - make pdebug -j 8 STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_LAPACK=${MFEM_USE_LAPACK:-"NO"} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" + make pdebug -j 8 STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_LAPACK=${MFEM_USE_LAPACK:-"NO"} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" MFEM_USE_SUPERLU=${MFEM_USE_SUPERLU:-"NO"} SUPERLU_DIR="$SUPERLU_DIR" SUPERLU_OPT="$SUPERLU_OPT" SUPERLU_LIB="$SUPERLU_LIB" check_result $? mfem-debug-installation fi cd $LIB_DIR @@ -132,7 +169,7 @@ else # NOTE(kevin): v4.5.2-dev commit. This is the mfem version used by PyMFEM v4.5.2.0. # This version matching is required to support pylibROM-PyMFEM interface. git checkout 00b2a0705f647e17a1d4ffcb289adca503f28d42 - make parallel -j 8 STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_LAPACK=${MFEM_USE_LAPACK:-"NO"} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" + make parallel -j 8 STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_LAPACK=${MFEM_USE_LAPACK:-"NO"} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" MFEM_USE_SUPERLU=${MFEM_USE_SUPERLU:-"NO"} SUPERLU_DIR="$SUPERLU_DIR" SUPERLU_OPT="$SUPERLU_OPT" SUPERLU_LIB="$SUPERLU_LIB" check_result $? mfem-parallel-installation fi cd $LIB_DIR From 6400aa50a8dc54ea66e2493da3eeb2e45c0630bb Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Mon, 25 Mar 2024 10:03:07 -0700 Subject: [PATCH 26/31] Fix getSpatialBasis calls to remove time parameter --- examples/prom/elliptic_eigenproblem_global_rom.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 33ebb13c8..4657a4dc0 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -479,11 +479,11 @@ int main(int argc, char *argv[]) const CAROM::Matrix *spatialbasis; if (rdim != -1) { - spatialbasis = reader.getSpatialBasis(0.0, rdim); + spatialbasis = reader.getSpatialBasis(rdim); } else { - spatialbasis = reader.getSpatialBasis(0.0, ef); + spatialbasis = reader.getSpatialBasis(ef); } const int numRowRB = spatialbasis->numRows(); From 86c74d1d8bfb1ac3840f0c7a55b628a13f20bd59 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Mon, 25 Mar 2024 12:27:58 -0700 Subject: [PATCH 27/31] Move helper functions to utilities --- .../prom/elliptic_eigenproblem_global_rom.cpp | 51 ++++--------------- lib/mfem/Utilities.cpp | 35 +++++++++++++ lib/mfem/Utilities.hpp | 44 ++++++++++++++++ 3 files changed, 89 insertions(+), 41 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 4657a4dc0..71e998529 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -833,43 +833,12 @@ double Potential(const Vector &x) const double min_d = 2.0 * sqrt(2.0); - auto check_within_portion = [&](const Vector &t, const double limit) { - // helper function to check if t is within limit percentage relative to the center of the mesh - for (int i = 0; i < dim; i++) - { - double domain_limit = limit * (bb_max[i] - bb_min[i]); - double mesh_center = 0.5 * (bb_max[i] + bb_min[i]); - - // check that t is within the limit relative to the center of the mesh - if ((t(i) - mesh_center) - (0.5 * domain_limit) > 1.0e-14) - { - std::cerr << "Error: value of t exceeds domain limit: t = " << t( - i) << ", limit = " << 0.5 * domain_limit << "\n"; - exit(-1); - } - } - }; - - auto map_to_ref_mesh = [](const double &bb_min, const double &bb_max, - const double &fraction) -> double { - // helper function to map a fractional value from [-1, 1] to [bb_min, bb_max] - CAROM_VERIFY(fraction <= 1.0 && fraction >= -1.0); - return bb_min + (fraction + 1.0) * ((bb_max - bb_min) * 0.5); - }; - - auto map_from_ref_mesh = [](const double &bb_min, const double &bb_max, - const double &value) -> double { - // helper function to map a value from the mesh range [bb_min, bb_max] to [-1, 1] - CAROM_VERIFY(value <= bb_max && value >= bb_min); - return -1.0 + (value - bb_min) * ((2.0) / (bb_max - bb_min)); - }; - X = x; for (int i = 0; i < dim; i++) { // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) - center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], center(i)); - neg_center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i)); + center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], center(i)); + neg_center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i)); } switch (problem) @@ -897,8 +866,8 @@ double Potential(const Vector &x) for (int i = 0; i < dim; i++) { // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) - center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], center(i)); - neg_center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i)); + center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], center(i)); + neg_center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i)); } return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp( -X.DistanceSquaredTo(neg_center) / c)); @@ -908,7 +877,7 @@ double Potential(const Vector &x) // min_d is the lower bound of the atom distance over time // verify t is within inner 20% of domain (centered around mesh origin) - check_within_portion(center, 0.2); + CAROM::check_within_portion(bb_min, bb_max, center, 0.2); return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d, 2)) + std::exp( @@ -924,8 +893,8 @@ double Potential(const Vector &x) center(1) += 0.5 * sin(2.0*kappa); // map the absolute location back to a fraction of the mesh domain - center(0) = map_from_ref_mesh(bb_min[0], bb_max[0], center(0)); - center(1) = map_from_ref_mesh(bb_min[1], bb_max[1], center(1)); + center(0) = CAROM::map_from_ref_mesh(bb_min[0], bb_max[0], center(0)); + center(1) = CAROM::map_from_ref_mesh(bb_min[1], bb_max[1], center(1)); neg_center = center; neg_center.Neg(); @@ -933,12 +902,12 @@ double Potential(const Vector &x) for (int i = 0; i < dim; i++) { // map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5) - center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], center(i)); - neg_center(i) = map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i)); + center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], center(i)); + neg_center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i)); } // verify t is within inner 20% of domain (centered around mesh origin) - check_within_portion(center, 0.2); + CAROM::check_within_portion(bb_min, bb_max, center, 0.2); return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d, 2)) + std::exp( diff --git a/lib/mfem/Utilities.cpp b/lib/mfem/Utilities.cpp index ebef43d84..6fcf4330b 100644 --- a/lib/mfem/Utilities.cpp +++ b/lib/mfem/Utilities.cpp @@ -62,4 +62,39 @@ void ComputeCtAB_vec(const HypreParMatrix& A, C.transposeMult(AB_carom, CtAB_vec); } +void check_within_portion(const mfem::Vector &bb_min, + const mfem::Vector &bb_max, const mfem::Vector &t, const double limit) +{ + // helper function to check if t is within limit percentage relative to the center of the mesh + for (int i = 0; i < t.Size(); i++) + { + double domain_limit = limit * (bb_max[i] - bb_min[i]); + double mesh_center = 0.5 * (bb_max[i] + bb_min[i]); + + // check that t is within the limit relative to the center of the mesh + if ((t(i) - mesh_center) - (0.5 * domain_limit) > 1.0e-14) + { + std::cerr << "Error: value of t exceeds domain limit: t = " << t( + i) << ", limit = " << 0.5 * domain_limit << "\n"; + exit(-1); + } + } +} + +double map_to_ref_mesh(const double &bb_min, const double &bb_max, + const double &fraction) +{ + // helper function to map a fractional value from [-1, 1] to [bb_min, bb_max] + CAROM_VERIFY(fraction <= 1.0 && fraction >= -1.0); + return bb_min + (fraction + 1.0) * ((bb_max - bb_min) * 0.5); +} + +double map_from_ref_mesh(const double &bb_min, const double &bb_max, + const double &value) +{ + // helper function to map a value from the mesh range [bb_min, bb_max] to [-1, 1] + CAROM_VERIFY(value <= bb_max && value >= bb_min); + return -1.0 + (value - bb_min) * ((2.0) / (bb_max - bb_min)); +} + } // end namespace CAROM diff --git a/lib/mfem/Utilities.hpp b/lib/mfem/Utilities.hpp index 464d1c93f..6caa503fe 100644 --- a/lib/mfem/Utilities.hpp +++ b/lib/mfem/Utilities.hpp @@ -55,6 +55,50 @@ void ComputeCtAB_vec(const HypreParMatrix& A, const CAROM::Matrix& C, CAROM::Vector& CtAB_vec); +/** + * @brief Helper function to ensure that @p t is within a given percentage of the domain relative to the center of the mesh + * + * @param bb_min Minimum corner of mesh bounding box. See mfem::Mesh::GetBoundingBox(). + * + * @param bb_max Maximum corner of mesh bounding box. See mfem::Mesh::GetBoundingBox(). + * + * @param t Point to check if within given @p limit percentage of mesh domain relative to mesh center. + * + * @param limit Fractional percentage (from [0, 1]) of mesh domain to check bounds of @p t. + */ +void check_within_portion(const mfem::Vector &bb_min, + const mfem::Vector &bb_max, const mfem::Vector &t, const double limit); + +/** + * @brief Maps a value from [-1, 1] to the corresponding mesh domain [bb_min, bb_max] + * + * @param bb_min Minimum value of domain + * + * @param bb_max Maximum value of domain + * + * @param fraction Value between [-1, 1] to map + * + * @returns @p fraction mapped to [bb_min, bb_max] + * @see map_from_ref_mesh + */ +double map_to_ref_mesh(const double &bb_min, const double &bb_max, + const double &fraction); + +/** + * @brief Maps a value within mesh domain [bb_min, bb_max] to the corresponding value between [-1, 1] + * + * @param bb_min Minimum value of domain + * + * @param bb_max Maximum value of domain + * + * @param value Value between [bb_min, bb_max] to map + * + * @returns @p value mapped to [-1, 1] + * @see map_to_ref_mesh + */ +double map_from_ref_mesh(const double &bb_min, const double &bb_max, + const double &value); + } // namespace CAROM #endif // MFEMUTILITIES_H From d12511cddc0d66e2d54c05b8dbceb230a2acf05d Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Tue, 26 Mar 2024 09:12:03 -0700 Subject: [PATCH 28/31] Fix issues with lapack/superlu setup --- scripts/setup.sh | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/scripts/setup.sh b/scripts/setup.sh index fb9374625..899f28b6e 100755 --- a/scripts/setup.sh +++ b/scripts/setup.sh @@ -105,9 +105,14 @@ METIS_DIR=$LIB_DIR/parmetis-4.0.3 METIS_OPT=-I${METIS_DIR}/metis/include METIS_LIB="-L${METIS_DIR}/build/lib/libparmetis -lparmetis -L${METIS_DIR}/build/lib/libmetis -lmetis" +if [ ${MFEM_USE_LAPACK} == "On" ]; then + # MFEM makefile needs YES/NO + MFEM_USE_LAPACK="YES" +fi + # Install distributed version of SuperLU cd $LIB_DIR -if [ -z ${INSTALL_SUPERLU} ] && [ ! -d "superlu_dist" ]; then +if [ ! -z ${INSTALL_SUPERLU} ] && [ ! -d "superlu_dist" ]; then wget https://github.com/xiaoyeli/superlu_dist/archive/refs/tags/v6.3.1.tar.gz tar -zxvf v6.3.1.tar.gz mv superlu_dist-6.3.1 superlu_dist @@ -117,6 +122,7 @@ if [ -z ${INSTALL_SUPERLU} ] && [ ! -d "superlu_dist" ]; then export PARMETIS_ROOT=${METIS_DIR} export PARMETIS_BUILD_DIR=${PARMETIS_ROOT}/build/lib/ cmake -DCMAKE_INSTALL_PREFIX=./ \ + -DCMAKE_INSTALL_LIBDIR=lib \ -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_C_COMPILER=mpicc \ -DCMAKE_CXX_COMPILER=mpicxx \ @@ -125,7 +131,7 @@ if [ -z ${INSTALL_SUPERLU} ] && [ ! -d "superlu_dist" ]; then -Denable_examples=OFF \ -DTPL_ENABLE_PARMETISLIB=ON \ -DTPL_PARMETIS_INCLUDE_DIRS="${PARMETIS_ROOT}/include;${PARMETIS_ROOT}/metis/include" \ - -DTPL_PARMETIS_LIBRARIES="${PARMETIS_BUILD_DIR}/libparmetis/libparmetis.a;${PARMETIS_BUILD_DIR}/libmetis/libmetis.a" \ + -DTPL_PARMETIS_LIBRARIES="${PARMETIS_BUILD_DIR}/libparmetis/libparmetis.so;${PARMETIS_BUILD_DIR}/libmetis/libmetis.so" \ -DTPL_ENABLE_COMBBLASLIB=OFF \ -DUSE_XSDK_DEFAULTS="False" \ -DBUILD_SHARED_LIBS=ON \ @@ -152,7 +158,7 @@ if [[ $BUILD_TYPE == "Debug" ]]; then if [[ $UPDATE_LIBS == "true" ]]; then cd mfem_debug git pull - make pdebug -j 8 STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_LAPACK=${MFEM_USE_LAPACK:-"NO"} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" MFEM_USE_SUPERLU=${MFEM_USE_SUPERLU:-"NO"} SUPERLU_DIR="$SUPERLU_DIR" SUPERLU_OPT="$SUPERLU_OPT" SUPERLU_LIB="$SUPERLU_LIB" + make -j 8 pdebug STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_LAPACK=${MFEM_USE_LAPACK:-"NO"} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" MFEM_USE_SUPERLU=${MFEM_USE_SUPERLU:-"NO"} SUPERLU_DIR="$SUPERLU_DIR" SUPERLU_OPT="$SUPERLU_OPT" SUPERLU_LIB="$SUPERLU_LIB" check_result $? mfem-debug-installation fi cd $LIB_DIR @@ -169,7 +175,7 @@ else # NOTE(kevin): v4.5.2-dev commit. This is the mfem version used by PyMFEM v4.5.2.0. # This version matching is required to support pylibROM-PyMFEM interface. git checkout 00b2a0705f647e17a1d4ffcb289adca503f28d42 - make parallel -j 8 STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_LAPACK=${MFEM_USE_LAPACK:-"NO"} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" MFEM_USE_SUPERLU=${MFEM_USE_SUPERLU:-"NO"} SUPERLU_DIR="$SUPERLU_DIR" SUPERLU_OPT="$SUPERLU_OPT" SUPERLU_LIB="$SUPERLU_LIB" + make -j 8 parallel STATIC=NO SHARED=YES MFEM_USE_MPI=YES MFEM_USE_GSLIB=${MG} MFEM_USE_LAPACK=${MFEM_USE_LAPACK:-"NO"} MFEM_USE_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_LIB" MFEM_USE_SUPERLU=${MFEM_USE_SUPERLU:-"NO"} SUPERLU_DIR="$SUPERLU_DIR" SUPERLU_OPT="$SUPERLU_OPT" SUPERLU_LIB="$SUPERLU_LIB" check_result $? mfem-parallel-installation fi cd $LIB_DIR From d0f7d25f2cffdd5bc826975ee936d30e3acc3bf8 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Tue, 26 Mar 2024 09:33:46 -0700 Subject: [PATCH 29/31] Add sample outputs for example problem --- .../prom/elliptic_eigenproblem_global_rom.cpp | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 71e998529..0167ec1a7 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -33,8 +33,30 @@ // // FOM run (for error calculation): // elliptic_eigenproblem_global_rom -fom -p 2 -f 1.0 +// Example output: +// Number of unknowns: 289 +// Eigenvalue 0: 0.04533314 +// Eigenvalue 1: 0.11332411 +// Eigenvalue 2: 0.11486387 +// Eigenvalue 3: 0.18192 +// Eigenvalue 4: 0.22964377 +// Elapsed time for assembling FOM: 1.471708e-03 second +// Elapsed time for solving FOM: 3.416310e-01 second // // Online phase: elliptic_eigenproblem_global_rom -online -p 2 -f 1.0 +// Example output: +// Eigenvalue 0: = 0.048430949 +// Eigenvalue 1: = 0.12021157 +// Eigenvalue 2: = 0.12147847 +// Eigenvalue 3: = 0.19456504 +// Eigenvalue 4: = 0.24285855 +// Relative error of ROM solution for eigenvalue 0 = 0.068334316 +// Relative error of ROM solution for eigenvalue 1 = 0.060776586 +// Relative error of ROM solution for eigenvalue 2 = 0.057586388 +// Relative error of ROM solution for eigenvalue 3 = 0.069508795 +// Relative error of ROM solution for eigenvalue 4 = 0.057544708 +// Elapsed time for assembling ROM: 4.289041e-03 second +// Elapsed time for solving ROM: 6.225410e-04 second #include "mfem.hpp" #include "linalg/BasisGenerator.h" From 92c3c4531a9826460b05be5219749446d38957f8 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Thu, 28 Mar 2024 10:01:25 -0700 Subject: [PATCH 30/31] Fix example sample commands --- examples/prom/elliptic_eigenproblem_global_rom.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index 0167ec1a7..b32b6e426 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -24,15 +24,15 @@ // order < 1 (quadratic for quadratic curvilinear mesh, NURBS for // NURBS mesh, etc.) // -// Offline phase: elliptic_eigenproblem_global_rom -offline -p 2 -id 0 -a 0.8 -// elliptic_eigenproblem_global_rom -offline -p 2 -id 1 -a 0.9 -// elliptic_eigenproblem_global_rom -offline -p 2 -id 2 -a 1.1 -// elliptic_eigenproblem_global_rom -offline -p 2 -id 3 -a 1.2 +// Offline phase: elliptic_eigenproblem_global_rom -offline -p 2 -id 0 -a 0.40 +// elliptic_eigenproblem_global_rom -offline -p 2 -id 1 -a 0.45 +// elliptic_eigenproblem_global_rom -offline -p 2 -id 2 -a 0.55 +// elliptic_eigenproblem_global_rom -offline -p 2 -id 3 -a 0.60 // // Merge phase: elliptic_eigenproblem_global_rom -merge -p 2 -ns 4 // // FOM run (for error calculation): -// elliptic_eigenproblem_global_rom -fom -p 2 -f 1.0 +// elliptic_eigenproblem_global_rom -fom -p 2 -a 0.5 // Example output: // Number of unknowns: 289 // Eigenvalue 0: 0.04533314 @@ -43,7 +43,7 @@ // Elapsed time for assembling FOM: 1.471708e-03 second // Elapsed time for solving FOM: 3.416310e-01 second // -// Online phase: elliptic_eigenproblem_global_rom -online -p 2 -f 1.0 +// Online phase: elliptic_eigenproblem_global_rom -online -p 2 -a 0.5 // Example output: // Eigenvalue 0: = 0.048430949 // Eigenvalue 1: = 0.12021157 From b8caa16b5a57b2e7e944654aad688d0531367f57 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Tue, 2 Apr 2024 10:36:28 -0700 Subject: [PATCH 31/31] Update function docs and fix abs error --- .../prom/elliptic_eigenproblem_global_rom.cpp | 4 +-- lib/mfem/Utilities.cpp | 12 ++++++--- lib/mfem/Utilities.hpp | 26 +++++++++++++------ 3 files changed, 28 insertions(+), 14 deletions(-) diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp index b32b6e426..a2087216c 100644 --- a/examples/prom/elliptic_eigenproblem_global_rom.cpp +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -899,7 +899,7 @@ double Potential(const Vector &x) // min_d is the lower bound of the atom distance over time // verify t is within inner 20% of domain (centered around mesh origin) - CAROM::check_within_portion(bb_min, bb_max, center, 0.2); + CAROM::verify_within_portion(bb_min, bb_max, center, 0.2); return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d, 2)) + std::exp( @@ -929,7 +929,7 @@ double Potential(const Vector &x) } // verify t is within inner 20% of domain (centered around mesh origin) - CAROM::check_within_portion(bb_min, bb_max, center, 0.2); + CAROM::verify_within_portion(bb_min, bb_max, center, 0.2); return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d, 2)) + std::exp( diff --git a/lib/mfem/Utilities.cpp b/lib/mfem/Utilities.cpp index 6fcf4330b..1c4ee424a 100644 --- a/lib/mfem/Utilities.cpp +++ b/lib/mfem/Utilities.cpp @@ -62,17 +62,21 @@ void ComputeCtAB_vec(const HypreParMatrix& A, C.transposeMult(AB_carom, CtAB_vec); } -void check_within_portion(const mfem::Vector &bb_min, - const mfem::Vector &bb_max, const mfem::Vector &t, const double limit) +void verify_within_portion(const mfem::Vector &bb_min, + const mfem::Vector &bb_max, + const mfem::Vector &t, const double limit) { - // helper function to check if t is within limit percentage relative to the center of the mesh + // helper function to check if t is within limit percentage relative + // to the center of the mesh + CAROM_VERIFY(t.Size() == bb_min.Size() && bb_min.Size() == bb_max.Size()); + CAROM_VERIFY(limit >= 0.0 && limit <= 1.0); for (int i = 0; i < t.Size(); i++) { double domain_limit = limit * (bb_max[i] - bb_min[i]); double mesh_center = 0.5 * (bb_max[i] + bb_min[i]); // check that t is within the limit relative to the center of the mesh - if ((t(i) - mesh_center) - (0.5 * domain_limit) > 1.0e-14) + if (std::abs((t(i) - mesh_center)) - (0.5 * domain_limit) > 1.0e-14) { std::cerr << "Error: value of t exceeds domain limit: t = " << t( i) << ", limit = " << 0.5 * domain_limit << "\n"; diff --git a/lib/mfem/Utilities.hpp b/lib/mfem/Utilities.hpp index 6caa503fe..9e7b32b57 100644 --- a/lib/mfem/Utilities.hpp +++ b/lib/mfem/Utilities.hpp @@ -56,18 +56,27 @@ void ComputeCtAB_vec(const HypreParMatrix& A, CAROM::Vector& CtAB_vec); /** - * @brief Helper function to ensure that @p t is within a given percentage of the domain relative to the center of the mesh + * @brief Helper function to ensure that @p t is within a given percentage of + * the domain relative to the center of the mesh. Performs the check for each + * dimension of the mesh (works if mesh is 2D or 3D). * - * @param bb_min Minimum corner of mesh bounding box. See mfem::Mesh::GetBoundingBox(). + * @param bb_min Minimum corner of mesh bounding box. + * See mfem::Mesh::GetBoundingBox(). * - * @param bb_max Maximum corner of mesh bounding box. See mfem::Mesh::GetBoundingBox(). + * @param bb_max Maximum corner of mesh bounding box. + * See mfem::Mesh::GetBoundingBox(). * - * @param t Point to check if within given @p limit percentage of mesh domain relative to mesh center. + * @param t Point to check if within given @p limit percentage of mesh domain + * relative to mesh center. * - * @param limit Fractional percentage (from [0, 1]) of mesh domain to check bounds of @p t. + * @param limit Fractional percentage (from [0, 1]) of mesh domain to check + * bounds of @p t. + * + * @note This will throw an error and exit if the check fails. */ -void check_within_portion(const mfem::Vector &bb_min, - const mfem::Vector &bb_max, const mfem::Vector &t, const double limit); +void verify_within_portion(const mfem::Vector &bb_min, + const mfem::Vector &bb_max, + const mfem::Vector &t, const double limit); /** * @brief Maps a value from [-1, 1] to the corresponding mesh domain [bb_min, bb_max] @@ -85,7 +94,8 @@ double map_to_ref_mesh(const double &bb_min, const double &bb_max, const double &fraction); /** - * @brief Maps a value within mesh domain [bb_min, bb_max] to the corresponding value between [-1, 1] + * @brief Maps a value within mesh domain [bb_min, bb_max] to the corresponding + * value between [-1, 1] * * @param bb_min Minimum value of domain *