From 336ae9801b850d10c0a3a0608e3de31fedbe48bc Mon Sep 17 00:00:00 2001 From: William Anderson <53445030+andersonw1@users.noreply.github.com> Date: Tue, 2 Apr 2024 11:17:43 -0700 Subject: [PATCH 1/2] Parametric dmdc heat conduction (#264) * Updated CMakeUpdated CMakeUpdated CMakeUpdated CMake * Algo files * Added file and headers * Changes to header files * DMDc online/predict phases compile and give output * Cleanup comments/print statements * Started adding interpolation for controls * Closer to interpolating controls * Controls are always interpolated now * Controls correctly interpolated * Remove ds.store * Cleanup * Changed logic in matrixinterpolator * Updated documentation at beginning * Changed documentation at start * Changed example at start and renamed visit output directories * astyle * documentation and formatting * Formatting on parametricDMDc * Added directory for outputs * Deleting pointers/making sure pointers are deleted * Deleting controls * Astyle * Copyright year and documentation in ParametricDMDc * formatting changes * astyle * DMDc.cpp formatting --- CMakeLists.txt | 4 +- .../dmd/parametric_dmdc_heat_conduction.cpp | 1074 +++++++++++++++++ lib/algo/DMDc.cpp | 12 +- lib/algo/DMDc.h | 39 +- lib/algo/ParametricDMDc.h | 196 +++ .../manifold_interp/MatrixInterpolator.cpp | 27 +- 6 files changed, 1341 insertions(+), 11 deletions(-) create mode 100644 examples/dmd/parametric_dmdc_heat_conduction.cpp create mode 100644 lib/algo/ParametricDMDc.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 4acea3b82..246cfcd3d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -153,7 +153,8 @@ if (ENABLE_EXAMPLES) local_tw_csv local_dw_csv parametric_tw_csv - parametric_dw_csv) + parametric_dw_csv + parametric_dmdc_heat_conduction) set(example_directories prom prom @@ -178,6 +179,7 @@ if (ENABLE_EXAMPLES) dmd dmd dmd + dmd dmd) list(LENGTH examples len1) diff --git a/examples/dmd/parametric_dmdc_heat_conduction.cpp b/examples/dmd/parametric_dmdc_heat_conduction.cpp new file mode 100644 index 000000000..9f7346322 --- /dev/null +++ b/examples/dmd/parametric_dmdc_heat_conduction.cpp @@ -0,0 +1,1074 @@ +/****************************************************************************** + * + * 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: Parametric DMDc for Heat_Conduction (adapted from ex16p.cpp) +// +// Compile with: make parametric_dmdc_heat_conduction +// +// ================================================================================= +// +// Sample runs and results for DMDc: +// +// Example changing parameter in state vector: +// mpirun -np 8 parametric_dmdc_heat_conduction -s 1 -offline --kappa 0.5 -rdim 6 +// mpirun -np 8 parametric_dmdc_heat_conduction -s 1 -offline --kappa 0.75 -rdim 6 +// mpirun -np 8 parametric_dmdc_heat_conduction -s 1 -offline --kappa 1.25 -rdim 6 +// mpirun -np 8 parametric_dmdc_heat_conduction -s 1 -offline --kappa 1.5 -rdim 6 +// +// mpirun -np 8 parametric_dmdc_heat_conduction -s 1 --kappa 1 -online -predict -visit +// +// Output: +// Relative error of DMDc temperature (u) at t_final: 0.5 is 0.0025944158 +// Elapsed time for solving FOM: 1.006058e+00 second +// Elapsed time for training DMDc: 7.535820e-01 second +// Elapsed time for predicting DMDc: 3.593333e-03 second +// +// Example changing parameter in control vector: +// mpirun -np 8 parametric_dmdc_heat_conduction -s 1 -offline -amp-in 1 -rdim 6 +// mpirun -np 8 parametric_dmdc_heat_conduction -s 1 -offline -amp-in 2 -rdim 6 +// mpirun -np 8 parametric_dmdc_heat_conduction -s 1 -offline -amp-in 4 -rdim 6 +// mpirun -np 8 parametric_dmdc_heat_conduction -s 1 -offline -amp-in 5 -rdim 6 +// +// mpirun -np 8 parametric_dmdc_heat_conduction -s 1 -amp-in 3 -online -predict -visit +// +// Output: +// Relative error of DMDc temperature (u) at t_final: 0.5 is 0.0028885833 +// Elapsed time for solving FOM: 9.007740e-01 second +// Elapsed time for training DMDc: 7.236116e-01 second +// Elapsed time for predicting DMDc: 5.489291e-03 second +// +// +// ================================================================================= +// +// Description: This example applies parametric DMDc to a time-dependent nonlinear +// heat equation problem of the form du/dt = C(u) + f. We use a +// nonlinear diffusion operator of the form +// C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u. We also assume +// time-varying external inlet and outlet sources which are +// respectively located at (0,0) and (0.5,0.5) in the reference +// domain [-1,1]^d. The sources have constant strengths given by +// amp-in and amp-out. + +// In this example, we can vary parameters in both the diffusion +// operator (kappa and alpha), and the control (amp-in and amp-out). +// In the offline stage, we first choose a set of training parameters +// and build individual DMDc's for each realization of the training +// parameters. After building the DMDc's, we then move to the online +// and predict phases. For a new set of parameter realizations, we +// build a new DMDc by interpolating the matrices Atilde, Btilde, and +// the controls from our training parameters. We then make predictions +// with our new DMDc and compare the speed and accuracy of our ROM +// with the full-order model. +// +// 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 "algo/DMDc.h" +#include "linalg/Vector.h" +#include +#include +#include +#include "algo/ParametricDMDc.h" +#include "utils/CSVDatabase.h" +#include "utils/HDFDatabase.h" + +#ifndef _WIN32 +#include // mkdir +#else +#include // _mkdir +#define mkdir(dir, mode) _mkdir(dir) +#endif + +using namespace std; +using namespace mfem; + +/** After spatial discretization, the conduction model can be written as: + * + * du/dt = M^{-1}(-Ku) + * + * where u is the vector representing the temperature, M is the mass matrix, + * and K is the diffusion operator with diffusivity depending on u: + * (\kappa + \alpha u). + * + * Class ConductionOperator represents the right-hand side of the above ODE. + */ +class ConductionOperator : public TimeDependentOperator +{ +protected: + ParFiniteElementSpace &fespace; + Array ess_tdof_list; // this list remains empty for pure Neumann b.c. + + ParBilinearForm *M; + ParBilinearForm *K; + ParLinearForm *B; + + HypreParMatrix Mmat; + HypreParMatrix Kmat; + HypreParMatrix *T; // T = M + dt K + double current_dt; + + CGSolver M_solver; // Krylov solver for inverting the mass matrix M + HypreSmoother M_prec; // Preconditioner for the mass matrix M + + CGSolver T_solver; // Implicit solver for T = M + dt K + HypreSmoother T_prec; // Preconditioner for the implicit solver + + FunctionCoefficient &src_func; // Source function coefficient + double alpha, kappa; // Nonlinear parameters + + mutable Vector z; // auxiliary vector + HypreParVector *b; // source vector + +public: + ConductionOperator(ParFiniteElementSpace &f, FunctionCoefficient &s, + double alpha, double kappa, const Vector &u); + + virtual void Mult(const Vector &u, Vector &du_dt) const; + /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k. + This is the only requirement for high-order SDIRK implicit integration.*/ + virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k); + + /// Update the diffusion BilinearForm K using the given true-dof vector `u`. + void SetSourceTime(const double t); + + /// Update the diffusion BilinearForm K using the given true-dof vector `u`. + void SetParameters(const Vector &u); + + virtual ~ConductionOperator(); +}; + +const CAROM::Matrix* +createControlMatrix(std::vector snapshots) +{ + CAROM_VERIFY(snapshots.size() > 0); + CAROM_VERIFY(snapshots[0]->dim() > 0); + for (int i = 0 ; i < snapshots.size() - 1; i++) + { + CAROM_VERIFY(snapshots[i]->dim() == snapshots[i + 1]->dim()); + CAROM_VERIFY(snapshots[i]->distributed() == snapshots[i + 1]->distributed()); + } + + CAROM::Matrix* snapshot_mat = new CAROM::Matrix(snapshots[0]->dim(), + snapshots.size(), + snapshots[0]->distributed()); + + for (int i = 0; i < snapshots[0]->dim(); i++) + { + for (int j = 0; j < snapshots.size(); j++) + { + snapshot_mat->item(i, j) = snapshots[j]->item(i); + } + } + + return snapshot_mat; +} + +double InitialTemperature(const Vector &x); +double TimeWindowFunction(const double t, const double t_begin, + const double t_end); +double Amplitude(const double t, const int index); +double SourceFunction(const Vector &x, double t); + +Vector bb_min, bb_max; // Mesh bounding box +double amp_in = 0.2; +double t_end_in = 0.1; +double amp_out = 0.1; +double t_end_out = 0.3; +double dt = 1.0e-2; + +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 = ""; + int ser_ref_levels = 2; + int par_ref_levels = 1; + int order = 2; + int ode_solver_type = 1; + double t_final = 0.5; + double alpha = 1.0e-2; + double kappa = 0.5; + int rdim = 6; + bool visualization = true; + bool visit = false; + int vis_steps = 5; + bool adios2 = false; + bool offline = false; + bool online = false; + + double closest_rbf_val = 0.9; + bool predict = false; + bool save_dofs = false; + bool csvFormat = true; + const char *basename = ""; + const char *temp_io_dir = "./outputs_parametric_dmdc_heat_conduction"; + std::string io_dir; + + 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(&ode_solver_type, "-s", "--ode-solver", + "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t" + "\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4."); + args.AddOption(&t_final, "-tf", "--t-final", + "Final time; start time is 0."); + args.AddOption(&dt, "-dt", "--time-step", + "Time step."); + args.AddOption(&alpha, "-a", "--alpha", + "Alpha coefficient."); + args.AddOption(&kappa, "-k", "--kappa", + "Kappa coefficient offset."); + args.AddOption(&_in, "-amp-in", "--amplitude-in", + "Amplitude of inlet source at (0,0)."); + args.AddOption(&_out, "-amp-out", "--amplitude-out", + "Amplitude of outlet source at (0.5,0.5)."); + args.AddOption(&t_end_in, "-t-end-in", "--t-end-in", + "End time of inlet source at (0,0)."); + args.AddOption(&t_end_out, "-t-end-out", "--t-end-out", + "End time of outlet source at (0.5,0.5)."); + 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(&adios2, "-adios2", "--adios2-streams", "-no-adios2", + "--no-adios2-streams", + "Save data using adios2 streams."); + args.AddOption(&rdim, "-rdim", "--rdim", + "Reduced dimension for DMDc."); + 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(&closest_rbf_val, "-crv", "--crv", + "DMD Closest RBF Value."); + args.AddOption(&predict, "-predict", "--predict", "-no-predict", "--no-predict", + "Enable or disable DMD prediction."); + args.AddOption(&save_dofs, "-save", "--save", "-no-save", "--no-save", + "Enable or disable MFEM DOF solution snapshot files)."); + args.AddOption(&csvFormat, "-csv", "--csv", "-hdf", "--hdf", + "Use CSV or HDF format for files output by -save option."); + args.AddOption(&basename, "-out", "--outputfile-name", + "Name of the sub-folder to dump files within the run directory."); + args.AddOption(&temp_io_dir, "-io", "--io-dir-name", + "Name of the sub-folder to load/dump input/output files within the current directory."); + + args.Parse(); + if (!args.Good()) + { + args.PrintUsage(cout); + MPI_Finalize(); + return 1; + } + + if (myid == 0) + { + args.PrintOptions(cout); + } + + io_dir = temp_io_dir; + mkdir(io_dir.c_str(), 0777); + + string outputPath = "."; + if (save_dofs) + { + outputPath = "run"; + if (string(basename) != "") { + outputPath += "/" + string(basename); + } + if (myid == 0) + { + const char path_delim = '/'; + string::size_type pos = 0; + do { + pos = outputPath.find(path_delim, pos+1); + string subdir = outputPath.substr(0, pos); + mkdir(subdir.c_str(), 0777); + } + while (pos != string::npos); + } + } + + // 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; + if (mesh_file == "") + { + mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL)); + } + else + { + mesh = new Mesh(mesh_file, 1, 1); + } + int dim = mesh->Dimension(); + + mesh->GetBoundingBox(bb_min, bb_max, max(order, 1)); + + // 4. Define the ODE solver used for time integration. Several implicit + // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as + // explicit Runge-Kutta methods are available. + ODESolver *ode_solver; + bool explicit_solver; + switch (ode_solver_type) + { + // Implicit L-stable methods + case 1: + ode_solver = new BackwardEulerSolver; + explicit_solver = false; + break; + case 2: + ode_solver = new SDIRK23Solver(2); + explicit_solver = false; + break; + case 3: + ode_solver = new SDIRK33Solver; + explicit_solver = false; + break; + // Explicit methods + case 11: + ode_solver = new ForwardEulerSolver; + explicit_solver = true; + break; + case 12: + ode_solver = new RK2Solver(0.5); + explicit_solver = true; + break; // midpoint method + case 13: + ode_solver = new RK3SSPSolver; + explicit_solver = true; + break; + case 14: + ode_solver = new RK4Solver; + explicit_solver = true; + break; + case 15: + ode_solver = new GeneralizedAlphaSolver(0.5); + explicit_solver = true; + break; + // Implicit A-stable methods (not L-stable) + case 22: + ode_solver = new ImplicitMidpointSolver; + explicit_solver = false; + break; + case 23: + ode_solver = new SDIRK23Solver; + explicit_solver = false; + break; + case 24: + ode_solver = new SDIRK34Solver; + explicit_solver = false; + break; + default: + cout << "Unknown ODE solver type: " << ode_solver_type << '\n'; + delete mesh; + return 3; + } + + // 5. 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(); + } + + // 6. 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(); + } + + // 7. Define the vector finite element space representing the current and the + // initial temperature, u_ref. + H1_FECollection fe_coll(order, dim); + ParFiniteElementSpace fespace(pmesh, &fe_coll); + + int fe_size = fespace.GlobalTrueVSize(); + if (myid == 0) + { + cout << "Number of temperature unknowns: " << fe_size << endl; + } + + ParGridFunction u_gf(&fespace); + + // 8. Set the initial conditions for u. All boundaries are considered + // natural. + FunctionCoefficient u_0(InitialTemperature); + u_gf.ProjectCoefficient(u_0); + Vector u; + u_gf.GetTrueDofs(u); + + // 9. Initialize the conduction operator and the VisIt visualization. + FunctionCoefficient s(SourceFunction); + ConductionOperator oper(fespace, s, alpha, kappa, u); + + u_gf.SetFromTrueDofs(u); + { + ostringstream mesh_name, sol_name; + mesh_name << io_dir << "/parametric_dmdc_heat_conduction_" << to_string( + alpha) << "_" << to_string(kappa) << "_" << to_string(amp_in) << "_" << + to_string(amp_out) << "-mesh." << setfill('0') << setw( + 6) << myid; + sol_name << io_dir << "/parametric_dmdc_heat_conduction_" << to_string( + alpha) << "_" << to_string(kappa) << "_" << to_string(amp_in) << "_" << + to_string(amp_out) << "-init." << setfill('0') << setw( + 6) << myid; + ofstream omesh(mesh_name.str().c_str()); + omesh.precision(precision); + pmesh->Print(omesh); + ofstream osol(sol_name.str().c_str()); + osol.precision(precision); + u_gf.Save(osol); + } + + VisItDataCollection visit_dc(io_dir + + "/parametric_dmdc_Heat_Conduction_FOM_" + + to_string(alpha) + "_" + to_string(kappa) + "_" + to_string( + amp_in) + "_" + to_string(amp_out), pmesh); + visit_dc.RegisterField("temperature", &u_gf); + if (visit) + { + visit_dc.SetCycle(0); + visit_dc.SetTime(0.0); + visit_dc.Save(); + } + + // Optionally output a BP (binary pack) file using ADIOS2. This can be + // visualized with the ParaView VTX reader. +#ifdef MFEM_USE_ADIOS2 + ADIOS2DataCollection* adios2_dc = NULL; + if (adios2) + { + std::string postfix(mesh_file); + postfix.erase(0, std::string("../data/").size() ); + postfix += "_o" + std::to_string(order); + postfix += "_solver" + std::to_string(ode_solver_type); + const std::string collection_name = "heat_conduction_dmdc-p-" + postfix + ".bp"; + + adios2_dc = new ADIOS2DataCollection(MPI_COMM_WORLD, collection_name, pmesh); + adios2_dc->SetParameter("SubStreams", std::to_string(num_procs/2) ); + adios2_dc->RegisterField("temperature", &u_gf); + adios2_dc->SetCycle(0); + adios2_dc->SetTime(0.0); + adios2_dc->Save(); + } +#endif + + 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 + { + sout.precision(precision); + sout << "solution\n" << *pmesh << u_gf; + sout << flush; + } + } + + StopWatch fom_timer, dmd_training_timer, dmd_prediction_timer; + + fom_timer.Start(); + + // 10. Perform time-integration (looping over the time iterations, ti, with a + // time-step dt). + ode_solver->Init(oper); + double t = 0.0; + vector ts; + double f[2]; + + CAROM::Vector* init = NULL; + + CAROM::Database *db = NULL; + if (csvFormat) + db = new CAROM::CSVDatabase(); + else + db = new CAROM::HDFDatabase(); + + vector snap_list; + + fom_timer.Stop(); + + + CAROM::DMDc* dmd_u = NULL; + f[0] = Amplitude(t, 0); + f[1] = Amplitude(t, 1); + + int dim_c = 2; // control dim + std::vector d_controls; // vector to store controls + CAROM::Vector* control = new CAROM::Vector(f, dim_c, false); + d_controls.push_back(control); + + if (offline) + { + dmd_training_timer.Start(); + + // 11. Create DMDc object and take initial sample. + u_gf.SetFromTrueDofs(u); + f[0] = Amplitude(t, 0); + f[1] = Amplitude(t, 1); + + std::cout << "t = " << t << std::endl; + + dmd_u = new CAROM::DMDc(u.Size(), 2, dt); + dmd_u->takeSample(u.GetData(), t, f, false); + + dmd_training_timer.Stop(); + } + + if (online) + { + u_gf.SetFromTrueDofs(u); + init = new CAROM::Vector(u.GetData(), u.Size(), true); + } + + if (save_dofs && myid == 0) + { + if (csvFormat) + { + mkdir((outputPath + "/step0").c_str(), 0777); + db->putDoubleArray(outputPath + "/step0/sol.csv", u.GetData(), u.Size()); + } + else + { + db->create(outputPath + "/dmd_0.hdf"); + db->putDoubleArray("step0sol", u.GetData(), u.Size()); + } + } + + ts.push_back(t); + snap_list.push_back(0); + + bool last_step = false; + for (int ti = 1; !last_step; ti++) + { + fom_timer.Start(); + + if (t + dt >= t_final - dt/2) + { + last_step = true; + } + + // Assuming Euler method is used + if (explicit_solver) + { + oper.SetSourceTime(t); + } + else + { + oper.SetSourceTime(t + dt); + } + ode_solver->Step(u, t, dt); + + fom_timer.Stop(); + + + f[0] = Amplitude(t, 0); + f[1] = Amplitude(t, 1); + if (offline) + { + dmd_training_timer.Start(); + + u_gf.SetFromTrueDofs(u); + dmd_u->takeSample(u.GetData(), t, f, last_step); + + if (myid == 0) + { + std::cout << "Taking snapshot at: t = " << t << std::endl; + } + + dmd_training_timer.Stop(); + } + + if (save_dofs && myid == 0) + { + if (csvFormat) + { + mkdir((outputPath + "/step" + to_string(ti)).c_str(), 0777); + db->putDoubleArray(outputPath + "/step" + to_string(ti) + "/sol.csv", + u.GetData(), u.Size()); + } + else + db->putDoubleArray("step" + to_string(ti) + "sol", + u.GetData(), u.Size()); + } + + ts.push_back(t); + snap_list.push_back(ti); + if (!last_step) + { + CAROM::Vector* control = new CAROM::Vector(f, dim_c, false); + d_controls.push_back(control); + } + + if (last_step || (ti % vis_steps) == 0) + { + if (myid == 0) + { + cout << "step " << ti << ", t = " << t << endl; + } + + u_gf.SetFromTrueDofs(u); + if (visualization) + { + sout << "parallel " << num_procs << " " << myid << "\n"; + sout << "solution\n" << *pmesh << u_gf << flush; + } + + if (visit) + { + visit_dc.SetCycle(ti); + visit_dc.SetTime(t); + visit_dc.Save(); + } + +#ifdef MFEM_USE_ADIOS2 + if (adios2) + { + adios2_dc->SetCycle(ti); + adios2_dc->SetTime(t); + adios2_dc->Save(); + } +#endif + } + oper.SetParameters(u); + } + +#ifdef MFEM_USE_ADIOS2 + if (adios2) + { + delete adios2_dc; + } +#endif + + if (save_dofs && myid == 0) + { + if (csvFormat) + { + db->putDoubleVector(outputPath + "/tval.csv", ts, ts.size()); + db->putInteger(outputPath + "/numsnap", snap_list.size()); + db->putIntegerArray(outputPath + "/snap_list.csv", snap_list.data(), + snap_list.size()); + } + else + { + db->putDoubleVector("tval", ts, ts.size()); + db->putInteger("numsnap", snap_list.size()); + db->putInteger("snap_bound_size", 0); + db->putIntegerArray("snap_list", snap_list.data(), + snap_list.size()); + } + } + + // 12. Save the final solution in parallel. This output can be viewed later + // using GLVis: "glvis -np -m heat_conduction_dmdc-mesh -g heat_conduction_dmdc-final". + { + ostringstream sol_name; + sol_name << io_dir << "/parametric_dmdc_Heat_Conduction_" << to_string( + alpha) << "_" << to_string(kappa) << "_" << to_string(amp_in) << "_" << + to_string(amp_out) << "-final." << setfill('0') << setw( + 6) << myid; + ofstream osol(sol_name.str().c_str()); + osol.precision(precision); + u_gf.Save(osol); + } + + // 13. Calculate the DMDc modes. + if (offline || online) + { + if (offline) + { + + if (myid == 0) + { + std::cout << "Creating DMDc with rdim: " << rdim << std::endl; + } + + dmd_u->train(rdim); + + dmd_training_timer.Stop(); + + dmd_u->save(io_dir + "/" + to_string(alpha) + "_" + to_string( + kappa) + "_" + to_string(amp_in) + "_" + to_string(amp_out)); + + if (myid == 0) + { + std::ofstream fout; + fout.open(io_dir + "/parameters.txt", std::ios::app); + fout << alpha << " " << kappa << " " << amp_in << " " << amp_out << std::endl; + fout.close(); + } + + const CAROM::Matrix* control_mat = createControlMatrix(d_controls); + + std::string full_file_name; + + full_file_name = io_dir + "/" + to_string(alpha) + "_" + to_string( + kappa) + "_" + to_string( + amp_in) + "_" + to_string(amp_out) + "_control"; + control_mat->write(full_file_name); + } + + + if (online) + { + if (myid == 0) + { + std::cout << "Creating DMD using the rdim of the offline phase" << std::endl; + } + + std::fstream fin(io_dir + "/parameters.txt", std::ios_base::in); + double curr_param; + std::vector dmdc_paths; + std::vector controls; + std::vector param_vectors; + + while (fin >> curr_param) + { + double curr_alpha = curr_param; + fin >> curr_param; + double curr_kappa = curr_param; + fin >> curr_param; + double curr_amp_in = curr_param; + fin >> curr_param; + double curr_amp_out = curr_param; + + dmdc_paths.push_back(io_dir + "/" + to_string(curr_alpha) + "_" + to_string( + curr_kappa) + "_" + to_string(curr_amp_in) + "_" + to_string(curr_amp_out) ); + + CAROM::Matrix* curr_control = new CAROM::Matrix(); + curr_control->read(io_dir + "/" + to_string(curr_alpha) + "_" + to_string( + curr_kappa) + "_" + to_string(curr_amp_in) + "_" + to_string( + curr_amp_out) + "_control"); + controls.push_back(curr_control); + + CAROM::Vector* param_vector = new CAROM::Vector(4, false); + param_vector->item(0) = curr_alpha; + param_vector->item(1) = curr_kappa; + param_vector->item(2) = curr_amp_in; + param_vector->item(3) = curr_amp_out; + param_vectors.push_back(param_vector); + } + fin.close(); + + CAROM::Vector* desired_param = new CAROM::Vector(4, false); + desired_param->item(0) = alpha; + desired_param->item(1) = kappa; + desired_param->item(2) = amp_in; + desired_param->item(3) = amp_out; + + dmd_training_timer.Start(); + + + CAROM::Matrix* controls_interpolated = new CAROM::Matrix(); + + CAROM::getParametricDMDc(dmd_u, param_vectors, dmdc_paths, controls, + controls_interpolated, desired_param, "G", "LS", closest_rbf_val); + + dmd_u->project(init,controls_interpolated); + + + dmd_training_timer.Stop(); + + delete desired_param; + delete controls_interpolated; + for (auto m : param_vectors) + delete m; + for (auto m : controls) + delete m; + } + + if (predict) + { + Vector true_solution_u(u.Size()); + true_solution_u = u.GetData(); + + dmd_prediction_timer.Start(); + + // 14. Predict the state at t_final using DMDc. + if (myid == 0) + { + std::cout << "Predicting temperature using DMDc" << std::endl; + } + + CAROM::Vector* result_u = dmd_u->predict(ts[0]); + Vector initial_dmd_solution_u(result_u->getData(), result_u->dim()); + u_gf.SetFromTrueDofs(initial_dmd_solution_u); + + VisItDataCollection dmd_visit_dc(io_dir + + "/parametric_dmdc_Heat_Conduction_ROM_" + + to_string(alpha) + "_" + to_string(kappa) + "_" + to_string( + amp_in) + "_" + to_string(amp_out), pmesh); + dmd_visit_dc.RegisterField("temperature", &u_gf); + + if (visit) + { + dmd_visit_dc.SetCycle(0); + dmd_visit_dc.SetTime(0.0); + dmd_visit_dc.Save(); + } + + delete result_u; + + if (visit) + { + for (int i = 1; i < ts.size(); i++) + { + if (i == ts.size() - 1 || (i % vis_steps) == 0) + { + result_u = dmd_u->predict(ts[i]); + Vector dmd_solution_u(result_u->getData(), result_u->dim()); + u_gf.SetFromTrueDofs(dmd_solution_u); + + dmd_visit_dc.SetCycle(i); + dmd_visit_dc.SetTime(ts[i]); + dmd_visit_dc.Save(); + + delete result_u; + } + } + } + + dmd_prediction_timer.Stop(); + + result_u = dmd_u->predict(t_final); + + // 15. Calculate the relative error between the DMDc final solution and the true solution. + Vector dmd_solution_u(result_u->getData(), result_u->dim()); + Vector diff_u(true_solution_u.Size()); + subtract(dmd_solution_u, true_solution_u, diff_u); + + double tot_diff_norm_u = sqrt(InnerProduct(MPI_COMM_WORLD, diff_u, diff_u)); + double tot_true_solution_u_norm = sqrt(InnerProduct(MPI_COMM_WORLD, + true_solution_u, true_solution_u)); + + if (myid == 0) + { + + std::cout << "Relative error of DMDc temperature (u) at t_final: " << t_final << + " is " << tot_diff_norm_u / tot_true_solution_u_norm << std::endl; + printf("Elapsed time for solving FOM: %e second\n", fom_timer.RealTime()); + printf("Elapsed time for training DMDc: %e second\n", + dmd_training_timer.RealTime()); + printf("Elapsed time for predicting DMDc: %e second\n", + dmd_prediction_timer.RealTime()); + } + + delete result_u; + + } + + } + + // 16. Free the used memory. + delete ode_solver; + delete pmesh; + delete dmd_u; + + MPI_Finalize(); + + return 0; +} + +ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, + FunctionCoefficient &s, + double al, double kap, const Vector &u) + : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), src_func(s), + M(NULL), K(NULL), T(NULL), current_dt(0.0), + M_solver(f.GetComm()), T_solver(f.GetComm()), z(height) +{ + const double rel_tol = 1e-8; + + M = new ParBilinearForm(&fespace); + M->AddDomainIntegrator(new MassIntegrator()); + M->Assemble(0); // keep sparsity pattern of M and K the same + M->FormSystemMatrix(ess_tdof_list, Mmat); + + M_solver.iterative_mode = false; + M_solver.SetRelTol(rel_tol); + M_solver.SetAbsTol(0.0); + M_solver.SetMaxIter(100); + M_solver.SetPrintLevel(0); + M_prec.SetType(HypreSmoother::Jacobi); + M_solver.SetPreconditioner(M_prec); + M_solver.SetOperator(Mmat); + + alpha = al; + kappa = kap; + + T_solver.iterative_mode = false; + T_solver.SetRelTol(rel_tol); + T_solver.SetAbsTol(0.0); + T_solver.SetMaxIter(100); + T_solver.SetPrintLevel(0); + T_solver.SetPreconditioner(T_prec); + + SetParameters(u); + + B = new ParLinearForm(&fespace); + B->AddDomainIntegrator(new DomainLFIntegrator(src_func)); +} + +void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const +{ + // Compute: + // du_dt = M^{-1}*-K(u) + // for du_dt + Kmat.Mult(u, z); + z.Neg(); // z = -z + z += *b; + M_solver.Mult(z, du_dt); +} + +void ConductionOperator::ImplicitSolve(const double dt, + const Vector &u, Vector &du_dt) +{ + // Solve the equation: + // du_dt = M^{-1}*[-K(u + dt*du_dt)] + // for du_dt + if (!T) + { + T = Add(1.0, Mmat, dt, Kmat); + current_dt = dt; + T_solver.SetOperator(*T); + } + MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt + Kmat.Mult(u, z); + z.Neg(); + z += *b; + T_solver.Mult(z, du_dt); +} + +void ConductionOperator::SetSourceTime(double t) +{ + src_func.SetTime(t); + B->Assemble(); + b = B->ParallelAssemble(); +} + +void ConductionOperator::SetParameters(const Vector &u) +{ + ParGridFunction u_alpha_gf(&fespace); + u_alpha_gf.SetFromTrueDofs(u); + for (int i = 0; i < u_alpha_gf.Size(); i++) + { + u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i); + } + + delete K; + K = new ParBilinearForm(&fespace); + + GridFunctionCoefficient u_coeff(&u_alpha_gf); + + K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff)); + K->Assemble(0); // keep sparsity pattern of M and K the same + K->FormSystemMatrix(ess_tdof_list, Kmat); + delete T; + T = NULL; // re-compute T on the next ImplicitSolve +} + +ConductionOperator::~ConductionOperator() +{ + delete T; + delete M; + delete K; + delete B; +} + +double InitialTemperature(const Vector &x) +{ + if (x.Norml2() < 0.5) + { + return 2.0; + } + else + { + return 1.0; + } +} + +double TimeWindowFunction(const double t, const double t_begin, + const double t_end) +{ + return 0.5 * (tanh((t - t_begin) / (5*dt)) - tanh((t - t_end) / (5*dt))); +} + +double Amplitude(const double t, const int index) +{ + if (index == 0) + { + return amp_in * TimeWindowFunction(t, 0.0, t_end_in); + } + else + { + return amp_out * TimeWindowFunction(t, 0.0, t_end_out); + } +} + +double SourceFunction(const Vector &x, double t) +{ + // map to the reference [-1,1] domain + Vector X(x.Size()), Y(x.Size()); + for (int i = 0; i < x.Size(); i++) + { + double center = (bb_min[i] + bb_max[i]) * 0.5; + X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]); + Y(i) = X(i) - 0.5; + } + + double r1 = X.Norml2() / 0.01; + double r2 = Y.Norml2() / 0.01; + return Amplitude(t,0) * exp(-0.5*r1*r1) - Amplitude(t,1) * exp(-0.5*r2*r2); +} + + diff --git a/lib/algo/DMDc.cpp b/lib/algo/DMDc.cpp index 1b1c0df53..7e4591767 100644 --- a/lib/algo/DMDc.cpp +++ b/lib/algo/DMDc.cpp @@ -107,7 +107,7 @@ DMDc::DMDc(std::string base_file_name) DMDc::DMDc(std::vector> eigs, Matrix* phi_real, Matrix* phi_imaginary, Matrix* B_tilde, int k, - double dt, double t_offset, Vector* state_offset) + double dt, double t_offset, Vector* state_offset, Matrix* basis) { // Get the rank of this process, and the number of processors. int mpi_init; @@ -128,6 +128,7 @@ DMDc::DMDc(std::vector> eigs, Matrix* phi_real, d_k = k; d_dt = dt; d_t_offset = t_offset; + d_basis = basis; setOffset(state_offset); } @@ -351,7 +352,7 @@ DMDc::constructDMDc(const Matrix* f_snapshots, d_sv.push_back(d_factorizer_in->S[i]); } - int d_k_in; + int d_k_in = d_k; if (d_energy_fraction != -1.0) { d_k_in = d_num_singular_vectors; @@ -404,6 +405,13 @@ DMDc::constructDMDc(const Matrix* f_snapshots, d_S_inv->item(i, i) = 1 / d_factorizer_in->S[static_cast(i)]; } + // Make sure the basis is freed since we are setting it with DMDc class instead + // of manually setting the basis when interpolating + if (d_basis) + { + delete d_basis; + } + if (B == NULL) { // SVD on outputs diff --git a/lib/algo/DMDc.h b/lib/algo/DMDc.h index e0360a8f2..f298c4cae 100644 --- a/lib/algo/DMDc.h +++ b/lib/algo/DMDc.h @@ -16,6 +16,7 @@ #ifndef included_DMDc_h #define included_DMDc_h +#include "ParametricDMDc.h" #include #include @@ -184,6 +185,41 @@ class DMDc void summary(std::string base_file_name); protected: + /** + * @brief Obtain DMD model interpolant at desired parameter point by + * interpolation of DMD models from training parameter points. + * + * @param[in] parametric_dmdc The interpolant DMD model at the desired point. + * @param[in] parameter_points The training parameter points. + * @param[in] dmdcs The DMD objects associated with + * each training parameter point. + * @param[in] controls The matrices of controls from previous + * runs which we use to interpolate. + * @param[in] controls_interpolated The interpolated controls. + * @param[in] desired_point The desired point at which to create a parametric DMD. + * @param[in] rbf The RBF type ("G" == gaussian, + * "IQ" == inverse quadratic, + * "IMQ" == inverse multiquadric) + * @param[in] interp_method The interpolation method type + * ("LS" == linear solve, + * "IDW" == inverse distance weighting, + * "LP" == lagrangian polynomials) + * @param[in] closest_rbf_val The RBF parameter determines the width of influence. + * Set the RBF value of the nearest two parameter points + * to a value between 0.0 to 1.0 + * @param[in] reorthogonalize_W Whether to reorthogonalize the interpolated W (basis) matrix. + */ + friend void getParametricDMDc(DMDc*& parametric_dmdc, + std::vector& parameter_points, + std::vector& dmdcs, + std::vector controls, + Matrix*& controls_interpolated, + Vector* desired_point, + std::string rbf, + std::string interp_method, + double closest_rbf_val, + bool reorthogonalize_W); + /** * @brief Constructor. Variant of DMDc with non-uniform time step size. * @@ -204,10 +240,11 @@ class DMDc * @param[in] dt d_dt * @param[in] t_offset d_t_offset * @param[in] state_offset d_state_offset + * @param[in] basis d_basis, set by DMDc class for offline stages. When interpolating a new DMDc, we enter the interpolated basis explicitly */ DMDc(std::vector> eigs, Matrix* phi_real, Matrix* phi_imaginary, Matrix* B_tilde, int k, - double dt, double t_offset, Vector* state_offset); + double dt, double t_offset, Vector* state_offset, Matrix* basis = nullptr); /** * @brief Unimplemented default constructor. diff --git a/lib/algo/ParametricDMDc.h b/lib/algo/ParametricDMDc.h new file mode 100644 index 000000000..78639a351 --- /dev/null +++ b/lib/algo/ParametricDMDc.h @@ -0,0 +1,196 @@ +/****************************************************************************** + * + * 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) + * + *****************************************************************************/ + +// Description: Computes the ParametricDMDc algorithm to obtain DMDc model interpolant +// at desired parameter point by interpolation of DMDc models at training parameter points. +// The implemented dynamic mode decomposition algorithm is derived from +// Tu et. al's paper "On Dynamic Mode Decomposition: Theory and +// Applications": https://arxiv.org/abs/1312.0041 +// The interpolation algorithm was adapted from "Gradient-based +// Constrained Optimization Using a Database of Linear Reduced-Order Models" +// by Y. Choi et al. We extend this algorithm to DMDc. + +#ifndef included_ParametricDMDc_h +#define included_ParametricDMDc_h + +#include "manifold_interp/MatrixInterpolator.h" +#include "linalg/Matrix.h" +#include "linalg/Vector.h" +#include "mpi.h" + +#include + +namespace CAROM { + +/** + * @brief Constructor. + * + * @param[in] parametric_dmdc The interpolant DMDc model at the desired point. + * @param[in] parameter_points The training parameter points. + * @param[in] dmdcs The DMDc objects associated with + * each training parameter point. + * @param[in] controls The matrices of controls from previous runs + * which we use to interpolate. + * @param[in] controls_interpolated The interpolated controls. + * @param[in] desired_point The desired point at which to create a parametric DMDc. + * @param[in] rbf The RBF type ("G" == gaussian, + * "IQ" == inverse quadratic, + * "IMQ" == inverse multiquadric) + * @param[in] interp_method The interpolation method type + * ("LS" == linear solve, + * "IDW" == inverse distance weighting, + * "LP" == lagrangian polynomials) + * @param[in] closest_rbf_val The RBF parameter determines the width of influence. + * Set the RBF value of the nearest two parameter points + * to a value between 0.0 to 1.0 + * @param[in] reorthogonalize_W Whether to reorthogonalize the interpolated W (basis) matrix. + */ +template +void getParametricDMDc(T*& parametric_dmdc, + std::vector& parameter_points, + std::vector& dmdcs, + std::vector controls, + Matrix*& controls_interpolated, + Vector* desired_point, + std::string rbf = "G", + std::string interp_method = "LS", + double closest_rbf_val = 0.9, + bool reorthogonalize_W = false) +{ + CAROM_VERIFY(parameter_points.size() == dmdcs.size()); + CAROM_VERIFY(dmdcs.size() > 1); + for (int i = 0; i < dmdcs.size() - 1; i++) + { + CAROM_VERIFY(dmdcs[i]->d_dt == dmdcs[i + 1]->d_dt); + CAROM_VERIFY(dmdcs[i]->d_k == dmdcs[i + 1]->d_k); + } + CAROM_VERIFY(closest_rbf_val >= 0.0 && closest_rbf_val <= 1.0); + + int mpi_init, rank; + MPI_Initialized(&mpi_init); + if (mpi_init == 0) { + MPI_Init(nullptr, nullptr); + } + + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + std::vector bases; + std::vector A_tildes; + std::vector B_tildes; + for (int i = 0; i < dmdcs.size(); i++) + { + bases.push_back(dmdcs[i]->d_basis); + A_tildes.push_back(dmdcs[i]->d_A_tilde); + B_tildes.push_back(dmdcs[i]->d_B_tilde); + } + + int ref_point = getClosestPoint(parameter_points, desired_point); + std::vector rotation_matrices = obtainRotationMatrices( + parameter_points, + bases, ref_point); + + CAROM::MatrixInterpolator basis_interpolator(parameter_points, + rotation_matrices, bases, ref_point, "B", rbf, interp_method, closest_rbf_val); + CAROM::Matrix* W = basis_interpolator.interpolate(desired_point, + reorthogonalize_W); + + CAROM::MatrixInterpolator A_tilde_interpolator(parameter_points, + rotation_matrices, A_tildes, ref_point, "R", rbf, interp_method, + closest_rbf_val); + CAROM::Matrix* A_tilde = A_tilde_interpolator.interpolate(desired_point); + + CAROM::MatrixInterpolator B_tilde_interpolator(parameter_points, + rotation_matrices, B_tildes, ref_point, "R", rbf, interp_method, + closest_rbf_val); + CAROM::Matrix* B_tilde = B_tilde_interpolator.interpolate(desired_point); + + CAROM::MatrixInterpolator control_interpolator(parameter_points, + rotation_matrices, controls, ref_point, "R", rbf, interp_method, + closest_rbf_val); + + controls_interpolated = control_interpolator.interpolate(desired_point); + + // Calculate the right eigenvalues/eigenvectors of A_tilde + ComplexEigenPair eigenpair = NonSymmetricRightEigenSolve(A_tilde); + std::vector> eigs = eigenpair.eigs; + + // Calculate phi (phi = W * eigenvectors) + CAROM::Matrix* phi_real = W->mult(eigenpair.ev_real); + CAROM::Matrix* phi_imaginary = W->mult(eigenpair.ev_imaginary); + + parametric_dmdc = new T(eigs, phi_real, phi_imaginary, B_tilde, + dmdcs[0]->d_k,dmdcs[0]->d_dt, + dmdcs[0]->d_t_offset, dmdcs[0]->d_state_offset, W); + + delete A_tilde; + delete eigenpair.ev_real; + delete eigenpair.ev_imaginary; + + for (auto m : rotation_matrices) + delete m; + +} + +/** + * @brief Constructor. + * + * @param[in] parameter_points The parameter points. + * @param[in] dmdc_paths The paths to the saved DMD objects associated with + * each parameter point. + * @param[in] desired_point The desired point at which to create a parametric DMD. + * @param[in] controls The matrices of controls from previous runs which we + * use to interpolate. + * @param[in] controls_interpolated The interpolated controls. + * @param[in] rbf The RBF type ("G" == gaussian, + * "IQ" == inverse quadratic, + * "IMQ" == inverse multiquadric) + * @param[in] interp_method The interpolation method type + * ("LS" == linear solve, + * "IDW" == inverse distance weighting, + * "LP" == lagrangian polynomials) + * @param[in] closest_rbf_val The RBF parameter determines the width of influence. + * Set the RBF value of the nearest two parameter points + * to a value between 0.0 to 1.0 + * @param[in] reorthogonalize_W Whether to reorthogonalize the interpolated W (basis) matrix. + */ +template +void getParametricDMDc(T*& parametric_dmdc, + std::vector& parameter_points, + std::vector& dmdc_paths, + std::vector controls, + Matrix*& controls_interpolated, + Vector* desired_point, + std::string rbf = "G", + std::string interp_method = "LS", + double closest_rbf_val = 0.9, + bool reorthogonalize_W = false) +{ + std::vector dmdcs; + for (int i = 0; i < dmdc_paths.size(); i++) + { + T* dmdc = new T(dmdc_paths[i]); + dmdcs.push_back(dmdc); + } + + getParametricDMDc(parametric_dmdc, parameter_points, dmdcs, controls, + controls_interpolated, + desired_point, + rbf, interp_method, closest_rbf_val, + reorthogonalize_W); + + for (int i = 0; i < dmdcs.size(); i++) + { + delete dmdcs[i]; + } +} + +} + +#endif diff --git a/lib/algo/manifold_interp/MatrixInterpolator.cpp b/lib/algo/manifold_interp/MatrixInterpolator.cpp index faa23b054..a0606344a 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.cpp +++ b/lib/algo/manifold_interp/MatrixInterpolator.cpp @@ -68,21 +68,34 @@ MatrixInterpolator::MatrixInterpolator(std::vector parameter_points, } else { - if (d_matrix_type == "B") + if (reduced_matrices[i]->numRows() == rotation_matrices[i]->numRows() + && reduced_matrices[i]->numColumns() == rotation_matrices[i]->numRows()) + { + Matrix* Q_tA = rotation_matrices[i]->transposeMult(reduced_matrices[i]); + Matrix* Q_tAQ = Q_tA->mult(rotation_matrices[i]); + delete Q_tA; + d_rotated_reduced_matrices.push_back(Q_tAQ); + d_rotated_reduced_matrices_owned.push_back(true); + } + else if (reduced_matrices[i]->numRows() == rotation_matrices[i]->numRows()) + { + Matrix* Q_tA = rotation_matrices[i]->transposeMult(reduced_matrices[i]); + d_rotated_reduced_matrices.push_back(Q_tA); + d_rotated_reduced_matrices_owned.push_back(true); + } + else if (reduced_matrices[i]->numColumns() == rotation_matrices[i]->numRows()) { Matrix* AQ = reduced_matrices[i]->mult(rotation_matrices[i]); d_rotated_reduced_matrices.push_back(AQ); + d_rotated_reduced_matrices_owned.push_back(true); } else { - Matrix* Q_tA = rotation_matrices[i]->transposeMult(reduced_matrices[i]); - Matrix* Q_tAQ = Q_tA->mult(rotation_matrices[i]); - delete Q_tA; - d_rotated_reduced_matrices.push_back(Q_tAQ); + d_rotated_reduced_matrices.push_back(reduced_matrices[i]); + d_rotated_reduced_matrices_owned.push_back(false); } - - d_rotated_reduced_matrices_owned.push_back(true); } + } } From e2b741674edcb35a4147da153379008362545a2c Mon Sep 17 00:00:00 2001 From: Coleman Kendrick <6309092+ckendrick@users.noreply.github.com> Date: Wed, 3 Apr 2024 16:53:34 -0600 Subject: [PATCH 2/2] Add elliptic eigenvalue global prom example (#266) * Add initial eigenvalue example * Rename example and update description * Add error calculation with FOM * Add options for energy fraction and reduced dim * Update eigenproblem problem setup * Fix center and kappa for problem 4 and 5 * Fix saving eigenvectors to visit * Rename example to elliptic eigenvalue * Update case 4 and 5 * Output initial condition and adjust center relative to mesh * Add case 6 and use square mesh by default * Update mesh bounds and mapping * Add additional case to vary center around point * Transform mesh to be centered around origin * Add modified cases for limiting domain * Refactor parameter mapping used in potential setup * Update case 9 setup to use absolute mesh location * Add eigenvector error calculation * Add timing info printout * Add option to compile mfem with lapack enabled * Update takeSample to use new function signature * Add verbose option for solver * Renaming mapping helper functions * Add SuperLU, Strumpack, MKL CPardiso preconditioners * Add SuperLU build to setup script * Fix getSpatialBasis calls to remove time parameter * Move helper functions to utilities * Fix issues with lapack/superlu setup * Add sample outputs for example problem * Fix example sample commands * Update function docs and fix abs error --- CMakeLists.txt | 2 + .../prom/elliptic_eigenproblem_global_rom.cpp | 939 ++++++++++++++++++ lib/mfem/Utilities.cpp | 39 + lib/mfem/Utilities.hpp | 54 + scripts/compile.sh | 12 +- scripts/setup.sh | 47 +- 6 files changed, 1090 insertions(+), 3 deletions(-) create mode 100644 examples/prom/elliptic_eigenproblem_global_rom.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 246cfcd3d..7a9d9ffa3 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 + elliptic_eigenproblem_global_rom dg_advection nonlinear_elasticity heat_conduction @@ -167,6 +168,7 @@ if (ENABLE_EXAMPLES) prom prom prom + prom dmd dmd dmd diff --git a/examples/prom/elliptic_eigenproblem_global_rom.cpp b/examples/prom/elliptic_eigenproblem_global_rom.cpp new file mode 100644 index 000000000..a2087216c --- /dev/null +++ b/examples/prom/elliptic_eigenproblem_global_rom.cpp @@ -0,0 +1,939 @@ +/****************************************************************************** + * + * 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: Elliptic Eigenproblem (adapted from ex11p.cpp) +// +// ================================================================================= +// +// 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: 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 -a 0.5 +// 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 -a 0.5 +// 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" +#include "linalg/BasisReader.h" +#include "linalg/Vector.h" +#include "mfem/Utilities.hpp" +#include +#include +#include + +using namespace std; +using namespace mfem; + +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[]) +{ + // 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 = ""; + 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; + double ef = 0.9999; + int rdim = -1; + int verbose_level = 0; + + int precision = 8; + cout.precision(precision); + + 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", + "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.AddOption(&ef, "-ef", "--energy_fraction", + "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."); +#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 (slu_solver && sp_solver) + { + 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) + { + 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; + 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, + 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 + { + 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(); + } + 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 + // 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. + assembleTimer.Start(); + ConstantCoefficient one(1.0); + + 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()) + { + ess_bdr.SetSize(pmesh->bdr_attributes.Max()); + ess_bdr = 1; + } + + ParBilinearForm *a = new ParBilinearForm(fespace); + a->AddDomainIntegrator(new DiffusionIntegrator(kappa_0)); + a->AddDomainIntegrator(new MassIntegrator(v_0)); + 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(); + +#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; + 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(verbose_level); + precond = amg; + } + else + { +#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); + lobpcg->SetPreconditioner(*precond); + lobpcg->SetMaxIter(200); + lobpcg->SetTol(1e-8); + lobpcg->SetPrecondUsageMode(1); + lobpcg->SetPrintLevel(verbose_level); + 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 + for (int i = 0; i < nev; i++) + { + if (myid == 0) + { + std::cout << " Eigenvalue " << i << ": " << eigenvalues[i] << "\n"; + } + if (offline) + { + x = lobpcg->GetEigenvector(i); + generator->takeSample(x.GetData()); + } + } + + if (offline) + { + 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; + if (rdim != -1) + { + spatialbasis = reader.getSpatialBasis(rdim); + } + else + { + spatialbasis = reader.getSpatialBasis(ef); + } + + 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 < nev; i++) + { + std::cout << "Eigenvalue " << i << ": = " << eigenvalues[i] << "\n"; + } + } + + 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; + 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; + } + + // 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 << "elliptic_eigenproblem-mesh." << setfill('0') << setw(6) << myid; + + ofstream mesh_ofs(mesh_name.str().c_str()); + 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) { + // convert eigenvector from HypreParVector to ParGridFunction + x = lobpcg->GetEigenvector(i); + } else { + // for online, eigenvectors are stored in evect matrix + Vector ev; + evect.GetRow(i, ev); + x = ev; + } + + 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(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(""); + } + + ofstream sol_ev_ofs(sol_ev_name.str().c_str()); + sol_ev_ofs.precision(16); + for (int i = 0; i < nev && i < eigenvalues.Size(); ++i) + { + sol_ev_ofs << eigenvalues[i] << std::endl; + } + } + + 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) + { + 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++) + { + if (fom || offline) { + // convert eigenvector from HypreParVector to ParGridFunction + x = lobpcg->GetEigenvector(i); + } else { + // for online, eigenvectors are stored in evect matrix + Vector ev; + evect.GetRow(i, ev); + x = ev; + } + 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; + 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 < nev && i < eigenvalues.Size(); i++) + { + if ( myid == 0 ) + { + cout << "Eigenmode " << i+1 << '/' << nev + << ", Lambda = " << eigenvalues[i] << endl; + } + + if (fom || offline) { + // convert eigenvector from HypreParVector to ParGridFunction + x = lobpcg->GetEigenvector(i); + } else { + // for online, eigenvectors are stored in evect matrix + Vector ev; + evect.GetRow(i, ev); + x = ev; + } + + sout << "parallel " << num_procs << " " << myid << "\n" + << "solution\n" << *pmesh << x << flush + << "window_title 'Eigenmode " << i+1 << '/' << nev + << ", Lambda = " << eigenvalues[i] << "'" << endl; + + char c; + if (myid == 0) + { + cout << "press (q)uit or (c)ontinue --> " << flush; + cin >> c; + } + MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD); + + if (c != 'c') + { + break; + } + } + sout.close(); + } + } + + // 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; + } + delete M; + delete A; + + delete fespace; + if (order > 0) + { + delete fec; + } + delete pmesh; + +#if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK) + delete Arow; +#endif + + MPI_Finalize(); + + return 0; +} + +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: + case 6: + case 7: + case 8: + case 9: + return 1.0; + } + return 0.0; +} + +double Potential(const Vector &x) +{ + int dim = x.Size(); + + 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(); + + // amplitude of gaussians for problems 4-6 + const double D = 100.0; + // width of gaussians for problems 4-6 + const double c = 0.05; + + const double min_d = 2.0 * sqrt(2.0); + + 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) = 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) + { + case 1: + case 2: + return 0.0; + case 3: + return 1.0; + case 4: + return D * std::exp(-X.DistanceSquaredTo(center) / c); + case 5: + return -D * std::exp(-X.DistanceSquaredTo(center) / c); + 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(); + + 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) = 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)); + 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) + 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( + -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 + + // 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) = 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(); + + 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) = 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) + 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( + -X.DistanceSquaredTo(neg_center) / std::pow(c * min_d, 2))); + } + return 0.0; +} diff --git a/lib/mfem/Utilities.cpp b/lib/mfem/Utilities.cpp index ebef43d84..1c4ee424a 100644 --- a/lib/mfem/Utilities.cpp +++ b/lib/mfem/Utilities.cpp @@ -62,4 +62,43 @@ void ComputeCtAB_vec(const HypreParMatrix& A, C.transposeMult(AB_carom, CtAB_vec); } +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 + 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 (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"; + 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..9e7b32b57 100644 --- a/lib/mfem/Utilities.hpp +++ b/lib/mfem/Utilities.hpp @@ -55,6 +55,60 @@ 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. 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_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. + * + * @note This will throw an error and exit if the check fails. + */ +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] + * + * @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 diff --git a/scripts/compile.sh b/scripts/compile.sh index 6dcc2243f..9e43a3a3f 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="Off" 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="On" + ;; 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} == "On" ]] && [[ ${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..899f28b6e 100755 --- a/scripts/setup.sh +++ b/scripts/setup.sh @@ -105,6 +105,49 @@ 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 + 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_INSTALL_LIBDIR=lib \ + -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.so;${PARMETIS_BUILD_DIR}/libmetis/libmetis.so" \ + -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 +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_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_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 @@ -132,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_METIS=YES MFEM_USE_METIS_5=YES METIS_DIR="$METIS_DIR" METIS_OPT="$METIS_OPT" METIS_LIB="$METIS_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