diff --git a/.gitignore b/.gitignore index 34dc5c0a6..016626da8 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,9 @@ dependencies *.log *.status +# Do not track backup files +*~ + # Do not track autogenerated config headers CAROM_config.h FCMangle.h diff --git a/examples/prom/mixed_nonlinear_diffusion.cpp b/examples/prom/mixed_nonlinear_diffusion.cpp index 2388d6cfc..8b6f7916f 100644 --- a/examples/prom/mixed_nonlinear_diffusion.cpp +++ b/examples/prom/mixed_nonlinear_diffusion.cpp @@ -31,16 +31,19 @@ // Analytic test (reproductive) // mpirun -n 1 ./mixed_nonlinear_diffusion -offline // mpirun -n 1 ./mixed_nonlinear_diffusion -merge -ns 1 -// mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -nsdim 20 -// mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -nsdim 20 -sopt +// mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -nsdim 20 -hrtype deim +// mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -nsdim 20 -hrtype qdeim +// mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -nsdim 20 -hrtype sopt // mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -nsdim 20 -ns 1 -eqp // // Relative l2 error of ROM solution at final timestep using DEIM sampling: 1.096776797994166e-08 -// Elapsed time for time integration loop using DEIM sampling: 0.6351594580000001 +// Elapsed time for time integration loop using DEIM sampling: 0.5686423310000001 +// Relative l2 error of ROM solution at final timestep using QDEIM sampling: 1.227055339859034e-08 +// Elapsed time for time integration loop using QDEIM sampling: 0.575228571 // Relative l2 error of ROM solution at final timestep using S_OPT sampling: 1.01945081122054e-08 -// Elapsed time for time integration loop using S_OPT sampling: 0.6669736559999999 -// Relative l2 error of ROM solution at final timestep using EQP: 1.46205848438194e-07 -// Elapsed time for time integration loop using EQP: 0.431521853 +// Elapsed time for time integration loop using S_OPT sampling: 0.5311964850000001 +// Relative l2 error of ROM solution at final timestep using EQP: 1.052559143494963e-08 +// Elapsed time for time integration loop using EQP: 0.346165296 // // Note that the timing of the time integration loop does not include setup, // which can be greater for S_OPT and EQP than for DEIM. @@ -48,16 +51,18 @@ // Initial step test (reproductive) // mpirun -n 1 ./mixed_nonlinear_diffusion -offline -p 1 // mpirun -n 1 ./mixed_nonlinear_diffusion -merge -ns 1 -p 1 -// mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -p 1 -// mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -p 1 -sopt +// mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -p 1 -hrtype deim +// mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -p 1 -hrtype sopt // mpirun -n 1 ./mixed_nonlinear_diffusion -online -rrdim 8 -rwdim 8 -nldim 20 -p 1 -ns 1 -eqp // // Relative l2 error of ROM solution at final timestep using DEIM sampling: 0.0003712362376412496 -// Elapsed time for time integration loop using DEIM sampling: 0.364855569 +// Elapsed time for time integration loop using DEIM sampling: 0.339288686 +// Relative l2 error of ROM solution at final timestep using QDEIM sampling: 0.000375066186385267 +// Elapsed time for time integration loop using QDEIM sampling: 0.32320409 // Relative l2 error of ROM solution at final timestep using S_OPT sampling: 0.0003797338657417907 -// Elapsed time for time integration loop using S_OPT sampling: 0.300462563 -// Relative l2 error of ROM solution at final timestep using EQP sampling: 0.0003710336208386964 -// Elapsed time for time integration loop using EQP sampling: 0.481740662 +// Elapsed time for time integration loop using S_OPT sampling: 0.32422266 +// Relative l2 error of ROM solution at final timestep using EQP sampling: 0.0003709977143117392 +// Elapsed time for time integration loop using EQP sampling: 0.457011311 // // Initial step parametric test (predictive) // mpirun -n 1 ./mixed_nonlinear_diffusion -p 1 -offline -id 0 -sh 0.25 @@ -65,16 +70,18 @@ // mpirun -n 1 ./mixed_nonlinear_diffusion -p 1 -offline -id 2 -sh 0.35 // mpirun -n 1 ./mixed_nonlinear_diffusion -p 1 -merge -ns 3 // mpirun -n 1 ./mixed_nonlinear_diffusion -p 1 -offline -id 3 -sh 0.3 -// mpirun -n 1 ./mixed_nonlinear_diffusion -p 1 -online -rrdim 8 -rwdim 8 -nldim 20 -sh 0.3 -id 3 -// mpirun -n 1 ./mixed_nonlinear_diffusion -p 1 -online -rrdim 8 -rwdim 8 -nldim 20 -sh 0.3 -id 3 -sopt +// mpirun -n 1 ./mixed_nonlinear_diffusion -p 1 -online -rrdim 8 -rwdim 8 -nldim 20 -sh 0.3 -id 3 -hrtype deim +// mpirun -n 1 ./mixed_nonlinear_diffusion -p 1 -online -rrdim 8 -rwdim 8 -nldim 20 -sh 0.3 -id 3 -hrtype sopt // mpirun -n 1 ./mixed_nonlinear_diffusion -p 1 -online -rrdim 8 -rwdim 8 -nldim 20 -sh 0.3 -id 3 -ns 3 -eqp -maxnnls 30 // // Relative l2 error of ROM solution at final timestep using DEIM sampling: 0.002681387312231006 -// Elapsed time for time integration loop using DEIM sampling: 0.355846074 +// Elapsed time for time integration loop using DEIM sampling: 0.359206509 +// Relative l2 error of ROM solution at final timestep using QDEIM sampling: 0.00266355670204476 +// Elapsed time for time integration loop using QDEIM sampling: 0.309432016 // Relative l2 error of ROM solution at final timestep using S_OPT sampling: 0.002701713369494112 -// Elapsed time for time integration loop using S_OPT sampling: 0.348985935 +// Elapsed time for time integration loop using S_OPT sampling: 0.41139788 // Relative l2 error of ROM solution at final timestep using EQP: 0.002659978000520714 -// Elapsed time for time integration loop using EQP sampling: 0.176821221 +// Elapsed time for time integration loop using EQP sampling: 0.164961644 // // Pointwise snapshots for DMD input // mpirun -n 1 ./mixed_nonlinear_diffusion -pwsnap -pwx 101 -pwy 101 @@ -90,9 +97,7 @@ #include "linalg/BasisGenerator.h" #include "linalg/BasisReader.h" #include "linalg/NNLS.h" -#include "hyperreduction/DEIM.h" -#include "hyperreduction/GNAT.h" -#include "hyperreduction/S_OPT.h" +#include "hyperreduction/Hyperreduction.h" #include "mfem/SampleMesh.hpp" #include "mfem/PointwiseSnapshot.hpp" @@ -550,7 +555,7 @@ int main(int argc, char *argv[]) bool offline = false; bool merge = false; bool online = false; - bool use_sopt = false; + const char *samplingType = "deim"; bool use_eqp = false; bool writeSampleMesh = false; int num_samples_req = -1; @@ -619,8 +624,8 @@ 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(&use_sopt, "-sopt", "--sopt", "-no-sopt", "--no-sopt", - "Use S-OPT sampling instead of DEIM for the hyperreduction."); + args.AddOption(&samplingType, "-hrtype", "--hrsamplingtype", + "Sampling type for hyperreduction."); args.AddOption(&num_samples_req, "-nsr", "--nsr", "Number of samples for the sampling algorithm to select."); args.AddOption(&use_eqp, "-eqp", "--eqp", "-no-eqp", "--no-eqp", @@ -950,8 +955,6 @@ int main(int argc, char *argv[]) CAROM::BasisReader readerFR("basisFR"); FR_librom = readerFR.getSpatialBasis(0.0); - // Compute sample points using DEIM, for hyperreduction - if (nldim == -1) { nldim = FR_librom->numColumns(); @@ -1005,6 +1008,7 @@ int main(int argc, char *argv[]) else { // Setup hyperreduction using DEIM, GNAT, or S-OPT + CAROM::Hyperreduction hr(samplingType); vector num_sample_dofs_per_proc(num_procs); if (num_samples_req != -1) @@ -1019,44 +1023,15 @@ int main(int argc, char *argv[]) // Now execute the chosen sampling algorithm to get the sampling information. Bsinv = new CAROM::Matrix(nsamp_R, nldim, false); vector sample_dofs(nsamp_R); // Indices of the sampled rows - if (use_sopt) - { - if (myid == 0) - printf("Using S_OPT sampling\n"); - CAROM::S_OPT(FR_librom, - nldim, - sample_dofs, - num_sample_dofs_per_proc, - *Bsinv, - myid, - num_procs, - nsamp_R); - } - else if (nsamp_R != nldim) - { - if (myid == 0) - printf("Using GNAT sampling\n"); - CAROM::GNAT(FR_librom, - nldim, - sample_dofs, - num_sample_dofs_per_proc, - *Bsinv, - myid, - num_procs, - nsamp_R); - } - else - { - if (myid == 0) - printf("Using DEIM sampling\n"); - CAROM::DEIM(FR_librom, - nldim, - sample_dofs, - num_sample_dofs_per_proc, - *Bsinv, - myid, - num_procs); - } + + hr.ComputeSamples(FR_librom, + nldim, + sample_dofs, + num_sample_dofs_per_proc, + *Bsinv, + myid, + num_procs, + nsamp_R); vector sample_dofs_S; // Indices of the sampled rows vector num_sample_dofs_per_proc_S(num_procs); @@ -1068,8 +1043,6 @@ int main(int argc, char *argv[]) readerS = new CAROM::BasisReader("basisS"); S_librom = readerS->getSpatialBasis(0.0); - // Compute sample points using DEIM - if (nsdim == -1) { nsdim = S_librom->numColumns(); @@ -1083,7 +1056,7 @@ int main(int argc, char *argv[]) if (myid == 0) printf("reduced S dim = %d\n",nsdim); - // Now execute the DEIM algorithm to get the sampling information. + // Now use a hyperreduction method to compute samples. if (num_samples_req != -1) { nsamp_S = num_samples_req; @@ -1095,38 +1068,16 @@ int main(int argc, char *argv[]) Ssinv = new CAROM::Matrix(nsamp_S, nsdim, false); sample_dofs_S.resize(nsamp_S); - if (use_sopt) - { - CAROM::S_OPT(S_librom, - nsdim, - sample_dofs_S, - num_sample_dofs_per_proc_S, - *Ssinv, - myid, - num_procs, - nsamp_S); - } - else if (nsamp_S != nsdim) - { - CAROM::GNAT(S_librom, - nsdim, - sample_dofs_S, - num_sample_dofs_per_proc_S, - *Ssinv, - myid, - num_procs, - nsamp_S); - } - else - { - CAROM::DEIM(S_librom, - nsdim, - sample_dofs_S, - num_sample_dofs_per_proc_S, - *Ssinv, - myid, - num_procs); - } + + + hr.ComputeSamples(S_librom, + nsdim, + sample_dofs_S, + num_sample_dofs_per_proc_S, + *Ssinv, + myid, + num_procs, + nsamp_S); } // Construct sample mesh diff --git a/examples/prom/nonlinear_elasticity_global_rom.cpp b/examples/prom/nonlinear_elasticity_global_rom.cpp index e004086a9..984f1624d 100644 --- a/examples/prom/nonlinear_elasticity_global_rom.cpp +++ b/examples/prom/nonlinear_elasticity_global_rom.cpp @@ -30,58 +30,72 @@ // and nonlinear term basis, with velocity initial condition: // // Offline phase: -// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 0.90 -id 0 +// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 0.9 -id 0 // -// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 1.10 -id 1 +// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 1.1 -id 1 // // Merge phase: // ./nonlinear_elasticity_global_rom -merge -ns 2 -dt 0.01 -tf 5.0 // // Create FOM comparison data: -// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 1.00 -id 2 +// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 1.0 -id 2 // // Online phase with full sampling: -// ./nonlinear_elasticity_global_rom -online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -rvdim 40 -rxdim 10 -hdim 71 -nsr 1170 -sc 1.00 +// ./nonlinear_elasticity_global_rom -online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -hrtype gnat -rvdim 40 -rxdim 10 -hdim 71 -nsr 1170 -sc 1.0 // Output message: // Elapsed time for time integration loop 1.80759 // Relative error of ROM position (x) at t_final: 5 is 0.000231698 // Relative error of ROM velocity (v) at t_final: 5 is 0.466941 // -// Online phase with strong hyper-reduction: -// ./nonlinear_elasticity_global_rom -online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -rvdim 40 -rxdim 10 -hdim 71 -nsr 100 -sc 1.00 +// Online phase with strong hyper-reduction, using GNAT (over-sampled DEIM): +// ./nonlinear_elasticity_global_rom -online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -hrtype gnat -rvdim 40 -rxdim 10 -hdim 71 -nsr 100 -sc 1.0 // Output message: -// Elapsed time for time integration loop 1.08048 +// Elapsed time for time integration loop 1.0111 // Relative error of ROM position (x) at t_final: 5 is 0.00209877 // Relative error of ROM velocity (v) at t_final: 5 is 1.39472 // +// Online phase with strong hyper-reduction, using QDEIM: +// ./nonlinear_elasticity_global_rom -online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -hrtype qdeim -rvdim 40 -rxdim 10 -hdim 71 -nsr 100 -sc 1.0 +// Output message: +// Elapsed time for time integration loop 1.02559 +// Relative error of ROM position (x) at t_final: 5 is 0.00188458 +// Relative error of ROM velocity (v) at t_final: 5 is 0.978726 +// // ================================================================================= // // Sample runs and results for parametric ROM using only displacement basis // and nonlinear term basis: // Offline phase: -// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 0.90 -xbo -def-ic -id 0 -// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 1.10 -xbo -def-ic -id 1 +// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 0.9 -xbo -def-ic -id 0 +// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 1.1 -xbo -def-ic -id 1 // // Merge phase: // ./nonlinear_elasticity_global_rom -merge -ns 2 -dt 0.01 -tf 5.0 -xbo // // Create FOM comparison data: -// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 1.00 -xbo -def-ic -id 2 +// ./nonlinear_elasticity_global_rom -offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 1.0 -xbo -def-ic -id 2 // // Online phase with full sampling: -// ./nonlinear_elasticity_global_rom -online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -rxdim 57 -hdim 183 -nsr 1170 -sc 1.00 -xbo -def-ic +// ./nonlinear_elasticity_global_rom -online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -hrtype gnat -rxdim 57 -hdim 183 -nsr 1170 -sc 1.0 -xbo -def-ic // Output message: // Elapsed time for time integration loop 18.9874 // Relative error of ROM position (x) at t_final: 5 is 7.08272e-05 // Relative error of ROM velocity (v) at t_final: 5 is 0.00387647 // -// Online phase with strong hyper reduction: -// ./nonlinear_elasticity_global_rom -online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -rxdim 2 -hdim 4 -nsr 10 -sc 1.00 -xbo -def-ic +// Online phase with strong hyper reduction, using GNAT (over-sampled DEIM): +// ./nonlinear_elasticity_global_rom -online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -hrtype gnat -rxdim 2 -hdim 4 -nsr 10 -sc 1.0 -xbo -def-ic // Output message: -// Elapsed time for time integration loop 1.01136 +// Elapsed time for time integration loop 0.120194 // Relative error of ROM position (x) at t_final: 5 is 0.0130818 // Relative error of ROM velocity (v) at t_final: 5 is 0.979978 // +// Online phase with strong hyper reduction, using QDEIM: +// ./nonlinear_elasticity_global_rom -online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -hrtype qdeim -rxdim 2 -hdim 4 -nsr 10 -sc 1.0 -xbo -def-ic +// Output message: +// Elapsed time for time integration loop 0.10806 +// Relative error of ROM position (x) at t_final: 5 is 0.0108709 +// Relative error of ROM velocity (v) at t_final: 5 is 1.30704 +// // This example runs in parallel with MPI, by using the same number of MPI ranks // in all phases (offline, merge, online). @@ -89,9 +103,7 @@ #include "linalg/Vector.h" #include "linalg/BasisGenerator.h" #include "linalg/BasisReader.h" -#include "hyperreduction/DEIM.h" -#include "hyperreduction/GNAT.h" -#include "hyperreduction/S_OPT.h" +#include "hyperreduction/Hyperreduction.h" #include "mfem/SampleMesh.hpp" #include "mfem/Utilities.hpp" @@ -166,7 +178,7 @@ class RomOperator : public TimeDependentOperator bool hyperreduce; bool x_base_only; - CAROM::Vector* pfom_librom, * pfom_v_librom; + CAROM::Vector *pfom_librom, *pfom_v_librom; Vector* pfom; Vector* pfom_x; Vector* pfom_v; @@ -262,7 +274,8 @@ CAROM::Matrix* GetFirstColumns(const int N, const CAROM::Matrix* A) // TODO: move this to the library? void BasisGeneratorFinalSummary(CAROM::BasisGenerator* bg, - const double energyFraction, int& cutoff, const std::string cutoffOutputPath) + const double energyFraction, int& cutoff, + const std::string cutoffOutputPath) { const int rom_dim = bg->getSpatialBasis()->numColumns(); const CAROM::Vector* sing_vals = bg->getSingularValues(); @@ -388,15 +401,15 @@ int main(int argc, char* argv[]) bool offline = false; bool merge = false; bool online = false; - bool use_sopt = false; bool hyperreduce = true; bool x_base_only = false; int num_samples_req = -1; + const char *samplingType = "gnat"; int nsets = 0; int id_param = 0; - // number of basis vectors to use + // Number of basis vectors to use int rxdim = -1; int rvdim = -1; int hdim = -1; @@ -444,8 +457,8 @@ 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(&use_sopt, "-sopt", "--sopt", "-no-sopt", "--no-sopt", - "Use S-OPT sampling instead of DEIM for the hyperreduction."); + args.AddOption(&samplingType, "-hrtype", "--hrsamplingtype", + "Sampling type for hyperreduction."); args.AddOption(&num_samples_req, "-nsr", "--nsr", "number of samples we want to select for the sampling algorithm."); args.AddOption(&rxdim, "-rxdim", "--rxdim", @@ -569,8 +582,8 @@ int main(int argc, char* argv[]) BlockVector vx(true_offset); ParGridFunction v_gf, x_gf; - v_gf.MakeTRef(&fespace, vx, - true_offset[0]); // Associate a new FiniteElementSpace and new true-dof data with the GridFunction. + // Associate a FiniteElementSpace and true-dof data with the GridFunctions. + v_gf.MakeTRef(&fespace, vx, true_offset[0]); x_gf.MakeTRef(&fespace, vx, true_offset[1]); ParGridFunction x_ref(&fespace); @@ -688,7 +701,7 @@ int main(int argc, char* argv[]) vis_v.precision(8); visualize(vis_v, pmesh, &x_gf, &v_gf, "Velocity", true); // Make sure all ranks have sent their 'v' solution before initiating - // another set of GLVis connections (one from each rank): + // another set of GLVis connections (one from each rank) MPI_Barrier(pmesh->GetComm()); vis_w.open(vishost, visport); if (vis_w) @@ -838,44 +851,17 @@ int main(int argc, char* argv[]) CAROM::Matrix* Hsinv = new CAROM::Matrix(nsamp_H, hdim, false); vector sample_dofs(nsamp_H); - if (use_sopt) - { - if (myid == 0) - printf("Using S_OPT sampling\n"); - CAROM::S_OPT(H_librom, - hdim, - sample_dofs, - num_sample_dofs_per_proc, - *Hsinv, - myid, - num_procs, - nsamp_H); - } - else if (nsamp_H != hdim) - { - if (myid == 0) - printf("Using GNAT sampling\n"); - CAROM::GNAT(H_librom, - hdim, - sample_dofs, - num_sample_dofs_per_proc, - *Hsinv, - myid, - num_procs, - nsamp_H); - } - else - { - if (myid == 0) - printf("Using DEIM sampling\n"); - CAROM::DEIM(H_librom, - hdim, - sample_dofs, - num_sample_dofs_per_proc, - *Hsinv, - myid, - num_procs); - } + + // Setup hyperreduction using DEIM, GNAT, or S-OPT + CAROM::Hyperreduction hr(samplingType); + hr.ComputeSamples(H_librom, + hdim, + sample_dofs, + num_sample_dofs_per_proc, + *Hsinv, + myid, + num_procs, + nsamp_H); // Construct sample mesh const int nspaces = 1; @@ -904,7 +890,8 @@ int main(int argc, char* argv[]) w_x = new CAROM::Vector(rxdim, false); *w = 0.0; - // Note that some of this could be done only on the ROM solver process, but it is tricky, since RomOperator assembles Bsp in parallel. + // Note that some of this could be done only on the ROM solver process, + // but it is tricky, since RomOperator assembles Bsp in parallel. wMFEM = new Vector(&((*w)(0)), rxdim + rvdim); // Initial conditions @@ -930,8 +917,8 @@ int main(int argc, char* argv[]) // 12. Set the initial conditions for v_gf, x_gf and vx, and define the // boundary conditions on a beam-like mesh (see description above). - sp_v_gf.MakeTRef(sp_XV_space, sp_vx, - sp_offset[0]); // Associate a new FiniteElementSpace and new true-dof data with the GridFunction. + // Associate a FiniteElementSpace and true-dof data with the GridFunctions. + sp_v_gf.MakeTRef(sp_XV_space, sp_vx, sp_offset[0]); sp_x_gf.MakeTRef(sp_XV_space, sp_vx, sp_offset[1]); VectorFunctionCoefficient* velo = 0; @@ -1117,16 +1104,16 @@ int main(int argc, char* argv[]) if (x_base_only == false && basis_generator_v->isNextSample(t)) { basis_generator_v->takeSample(vx_diff.GetBlock(0), t, dt); - basis_generator_v->computeNextSampleTime(vx_diff.GetBlock(0), dvdt.GetData(), - t); + basis_generator_v->computeNextSampleTime(vx_diff.GetBlock(0), + dvdt.GetData(), t); basis_generator_H->takeSample(oper.H_sp.GetData(), t, dt); } if (basis_generator_x->isNextSample(t)) { basis_generator_x->takeSample(vx_diff.GetBlock(1), t, dt); - basis_generator_x->computeNextSampleTime(vx_diff.GetBlock(1), dxdt.GetData(), - t); + basis_generator_x->computeNextSampleTime(vx_diff.GetBlock(1), + dxdt.GetData(), t); if (x_base_only == true) { @@ -1207,7 +1194,8 @@ int main(int argc, char* argv[]) if (offline) { - // Sample final solution, to prevent extrapolation in ROM between the last sample and the end of the simulation. + // Sample final solution, to prevent extrapolation in ROM between the + // last sample and the end of the simulation. dvxdt = oper.dvxdt_sp.GetData(); vx_diff = BlockVector(vx); vx_diff -= vx0; @@ -1303,10 +1291,10 @@ int main(int argc, char* argv[]) if (myid == 0) { - cout << "Relative error of ROM position (x) at t_final: " << t_final << - " is " << tot_diff_norm_x / tot_x_fom_norm << endl; - cout << "Relative error of ROM velocity (v) at t_final: " << t_final << - " is " << tot_diff_norm_v / tot_v_fom_norm << endl; + cout << "Relative error of ROM position (x) at t_final: " << t_final + << " is " << tot_diff_norm_x / tot_x_fom_norm << endl; + cout << "Relative error of ROM velocity (v) at t_final: " << t_final + << " is " << tot_diff_norm_v / tot_v_fom_norm << endl; } } diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 677d1a799..a523dd67d 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -50,6 +50,7 @@ set(module_list hyperreduction/S_OPT hyperreduction/STSampling hyperreduction/Utilities + hyperreduction/Hyperreduction utils/Database utils/HDFDatabase utils/CSVDatabase diff --git a/lib/hyperreduction/GNAT.cpp b/lib/hyperreduction/GNAT.cpp index 83c47b8f5..faa65a243 100644 --- a/lib/hyperreduction/GNAT.cpp +++ b/lib/hyperreduction/GNAT.cpp @@ -97,8 +97,8 @@ void GNAT(const Matrix* f_basis, const int ns_mod_nr = num_samples % num_basis_vectors; int ns = 0; - // The small matrix inverted by the algorithm. We'll allocate the largest - // matrix we'll need and set its size at each step in the algorithm. + // The small matrix inverted by the algorithm. We allocate the largest + // size needed and set the size at each step in the algorithm. Matrix M(num_samples, std::max(num_basis_vectors-1, 1), false); // Scratch space used throughout the algorithm. @@ -125,8 +125,8 @@ void GNAT(const Matrix* f_basis, std::vector all_num_init_samples(num_procs); std::vector all_init_samples(total_num_init_samples); - MPI_Allgather(&num_init_samples, 1, MPI_INT, all_num_init_samples.data(), 1, - MPI_INT, MPI_COMM_WORLD); + MPI_Allgather(&num_init_samples, 1, MPI_INT, all_num_init_samples.data(), + 1, MPI_INT, MPI_COMM_WORLD); for (int i = 0; i < myid; ++i) { @@ -134,7 +134,7 @@ void GNAT(const Matrix* f_basis, } } - // Figure out the 1st sampled rows of the RHS. + // Figure out the first sampled rows of the RHS. RowInfo f_bv_max_local, f_bv_max_global; const int ns0 = 0 < ns_mod_nr ? (num_samples / num_basis_vectors) + 1 : @@ -215,10 +215,14 @@ void GNAT(const Matrix* f_basis, double tmp = 0.0; for (int minv_col = 0; minv_col < ns; ++minv_col) { if (ns == i) + { tmp += M.item(minv_row, minv_col)*tmp_fs.item(minv_col, i); + } else - tmp += M.item(minv_col, minv_row)*tmp_fs.item(minv_col, - i); // Transposing M^+, which is stored as its transpose. + { + // Transposing M^+, which is stored as its transpose. + tmp += M.item(minv_col, minv_row)*tmp_fs.item(minv_col, i); + } } c[minv_row] = tmp; } @@ -288,7 +292,7 @@ void GNAT(const Matrix* f_basis, CAROM_ASSERT(num_samples == ns); - // Fill f_sampled_row, and f_sampled_rows_per_proc. Unscramble tmp_fs into + // Fill f_sampled_row, and f_sampled_rows_per_proc. Unscramble tmp_fs into // f_basis_sampled_inv. int idx = 0; for (int i = 0; i < num_procs; ++i) { diff --git a/lib/hyperreduction/Hyperreduction.cpp b/lib/hyperreduction/Hyperreduction.cpp new file mode 100644 index 000000000..b110dc326 --- /dev/null +++ b/lib/hyperreduction/Hyperreduction.cpp @@ -0,0 +1,87 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2023, 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) + * + *****************************************************************************/ + +#include "Hyperreduction.h" +#include "DEIM.h" +#include "QDEIM.h" +#include "GNAT.h" +#include "S_OPT.h" + +#include "linalg/Matrix.h" +#include "utils/Utilities.h" + +namespace CAROM { + +Hyperreduction::Hyperreduction(const char* sampling_type) +{ + auto iter = SamplingTypeMap.find(sampling_type); + CAROM_VERIFY(iter != std::end(SamplingTypeMap)); + + samplingType = iter->second; +} + +void Hyperreduction::ComputeSamples(const Matrix* f_basis, + int num_f_basis_vectors_used, + std::vector& f_sampled_row, + std::vector& f_sampled_rows_per_proc, + Matrix& f_basis_sampled_inv, + int myid, + int num_procs, + const int num_samples_req, + std::vector *init_samples, + bool qr_factorize) +{ + switch (samplingType) + { + case deim: + CAROM_VERIFY(num_samples_req == f_basis->numColumns()); + DEIM(f_basis, + num_f_basis_vectors_used, + f_sampled_row, + f_sampled_rows_per_proc, + f_basis_sampled_inv, + myid, num_procs); + return; + case gnat: + GNAT(f_basis, + num_f_basis_vectors_used, + f_sampled_row, + f_sampled_rows_per_proc, + f_basis_sampled_inv, + myid, num_procs, + num_samples_req, + init_samples); + return; + case qdeim: + QDEIM(f_basis, + num_f_basis_vectors_used, + f_sampled_row, + f_sampled_rows_per_proc, + f_basis_sampled_inv, + myid, num_procs, + num_samples_req); + return; + case sopt: + S_OPT(f_basis, + num_f_basis_vectors_used, + f_sampled_row, + f_sampled_rows_per_proc, + f_basis_sampled_inv, + myid, num_procs, + num_samples_req, + init_samples, + qr_factorize); + return; + default: + CAROM_ERROR("Sampling type not supported"); + } +} + +} diff --git a/lib/hyperreduction/Hyperreduction.h b/lib/hyperreduction/Hyperreduction.h new file mode 100644 index 000000000..4fd4fc91a --- /dev/null +++ b/lib/hyperreduction/Hyperreduction.h @@ -0,0 +1,73 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2023, 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: Interface to hyperreduction algorithms. + +#ifndef included_Hyperreduction_h +#define included_Hyperreduction_h + + +#include +#include + +#include + +namespace CAROM { + +enum SamplingType +{ + deim, // Default, DEIM + gnat, // GNAT + qdeim, // QDEIM + sopt // S-OPT +}; + +static std::unordered_map SamplingTypeMap = +{ + {"deim", deim}, + {"gnat", gnat}, + {"qdeim", qdeim}, + {"sopt", sopt} +}; + +class Matrix; + +class Hyperreduction +{ +public: + + Hyperreduction(SamplingType stype) : + samplingType(stype) + { } + + Hyperreduction(const char* sampling_type); + + void SetSamplingType(SamplingType stype) + { + samplingType = stype; + } + + void ComputeSamples(const Matrix* f_basis, + int num_f_basis_vectors_used, + std::vector& f_sampled_row, + std::vector& f_sampled_rows_per_proc, + Matrix& f_basis_sampled_inv, + int myid, + int num_procs, + const int num_samples_req = -1, + std::vector *init_samples=NULL, + bool qr_factorize = false); + +private: + SamplingType samplingType; +}; + +} +#endif diff --git a/lib/hyperreduction/QDEIM.cpp b/lib/hyperreduction/QDEIM.cpp index a07e66c65..34fbf3619 100644 --- a/lib/hyperreduction/QDEIM.cpp +++ b/lib/hyperreduction/QDEIM.cpp @@ -37,8 +37,8 @@ QDEIM(const Matrix* f_basis, // processor that owns each sampled row, and fills f_basis_sampled_inv with the // sampled rows of the basis of the RHS. - CAROM_VERIFY(num_f_basis_vectors_used == - f_basis->numColumns()); // The QR implementation uses the entire matrix. + // The QR implementation uses the entire matrix. + CAROM_VERIFY(num_f_basis_vectors_used == f_basis->numColumns()); CAROM_VERIFY(f_basis->numColumns() <= num_samples_req && num_samples_req <= f_basis->numDistributedRows()); CAROM_VERIFY(num_samples_req == f_basis_sampled_inv.numRows() @@ -99,7 +99,8 @@ QDEIM(const Matrix* f_basis, ns[owner]++; } - // Reorder f_sampled_row and f_sampled_row_owner to match the order of f_basis_sampled_inv + // Reorder f_sampled_row and f_sampled_row_owner to match the order + // of f_basis_sampled_inv for (int i=0; imult(V); // distributed CAROM_VERIFY(Ubt->distributed() && Ubt->numRows() == f_basis->numRows() @@ -254,9 +256,9 @@ QDEIM(const Matrix* f_basis, for (int i=0; iitem(i, j) * Ubt->item(i, - j); // column sums of Ub.^2, which are row sums of Ubt.^2 + r[i] += Ubt->item(i, j) * Ubt->item(i,j); r[i] -= sqrt((r[i] * r[i]) - (4.0 * g * Ubt->item(i, n-1) * Ubt->item(i, n-1))); } @@ -311,8 +313,9 @@ QDEIM(const Matrix* f_basis, globalSamples.insert(f_sampled_row[s]); } - // Send one row of f_basis, corresponding to sample s, to root process for f_basis_sampled_inv. - // First, scatter from root to tell the owning process the sample index. + // Send one row of f_basis, corresponding to sample s, to the root + // process for f_basis_sampled_inv. First, scatter from root to tell + // the owning process the sample index. int sample = -1; MPI_Scatter(ns.data(), 1, MPI_INT, &sample, 1, MPI_INT, 0, MPI_COMM_WORLD); @@ -321,22 +324,22 @@ QDEIM(const Matrix* f_basis, if (sample > -1) { CAROM_VERIFY(sample >= row_offset[myid]); - MPI_Send(f_basis->getData() + ((sample - row_offset[myid]) * numCol), numCol, - MPI_DOUBLE, 0, tagSendRecv, MPI_COMM_WORLD); + MPI_Send(f_basis->getData() + ((sample - row_offset[myid]) * numCol), + numCol, MPI_DOUBLE, 0, tagSendRecv, MPI_COMM_WORLD); } if (myid == 0) { MPI_Status status; - MPI_Recv(sampled_row_data.data() + (s*numCol), numCol, MPI_DOUBLE, owner, - tagSendRecv, MPI_COMM_WORLD, &status); + MPI_Recv(sampled_row_data.data() + (s*numCol), numCol, MPI_DOUBLE, + owner, tagSendRecv, MPI_COMM_WORLD, &status); } delete Ubt; } // loop s over samples - // Subtract row_offset to convert f_sampled_row from global to local indices - // Also, reorder f_sampled_row by process + // Subtract row_offset to convert f_sampled_row from global to local indices. + // Also, reorder f_sampled_row by process. if (myid == 0) { ns[0] = 0; @@ -414,8 +417,8 @@ QDEIM(const Matrix* f_basis, if (!f_basis->distributed()) { - CAROM_VERIFY(numCol == - num_samples_req); // GappyPOD+E not implemented for oversampling if not distributed + // GappyPOD+E not implemented for oversampling if not distributed + CAROM_VERIFY(numCol == num_samples_req); } // Now invert f_basis_sampled_inv, storing its transpose. diff --git a/lib/hyperreduction/S_OPT.cpp b/lib/hyperreduction/S_OPT.cpp index 7efa656b6..f21534dc0 100644 --- a/lib/hyperreduction/S_OPT.cpp +++ b/lib/hyperreduction/S_OPT.cpp @@ -380,8 +380,8 @@ S_OPT(const Matrix* f_basis, { tmp += g1.item(j, k) * GG.item(j, k); } - A->item(j) = std::max(0.0, ata + (Vo->item(j, i - 1) * Vo->item(j, - i - 1)) - tmp); + A->item(j) = std::max(0.0, ata + + (Vo->item(j, i - 1) * Vo->item(j, i - 1)) - tmp); } nV.setSize(i); @@ -491,7 +491,7 @@ S_OPT(const Matrix* f_basis, delete noM; } - // Fill f_sampled_row, and f_sampled_rows_per_proc. Unscramble tmp_fs into + // Fill f_sampled_row, and f_sampled_rows_per_proc. Unscramble tmp_fs into // f_basis_sampled_inv. int idx = 0; for (int i = 0; i < num_procs; ++i) { diff --git a/lib/hyperreduction/Utilities.h b/lib/hyperreduction/Utilities.h index 2fa04304d..2572a47c5 100644 --- a/lib/hyperreduction/Utilities.h +++ b/lib/hyperreduction/Utilities.h @@ -17,8 +17,6 @@ namespace CAROM { -#ifndef DOXYGEN_IGNORE - /** * @brief Struct to hold the local maximum absolute value of a basis vector, * the row it is in, and the processor that owns it. We will reduce this @@ -37,8 +35,6 @@ typedef struct */ void RowInfoMax(RowInfo* a, RowInfo* b, int* len, MPI_Datatype* type); -#endif /* DOXYGEN_IGNORE */ - } #endif diff --git a/regression_tests/tests/prom/mixed_nonlinear_diffusion.sh b/regression_tests/tests/prom/mixed_nonlinear_diffusion.sh index bf8af6063..035ac6f32 100755 --- a/regression_tests/tests/prom/mixed_nonlinear_diffusion.sh +++ b/regression_tests/tests/prom/mixed_nonlinear_diffusion.sh @@ -3,8 +3,9 @@ source $GITHUB_WORKSPACE/regression_tests/common.sh CMDS=( "./mixed_nonlinear_diffusion -offline -tf 0.01" "./mixed_nonlinear_diffusion -merge -ns 1 -tf 0.01" - "./mixed_nonlinear_diffusion -online -tf 0.01 -rrdim 2 -rwdim 2 -nldim 3" - "./mixed_nonlinear_diffusion -online -tf 0.01 -rrdim 2 -rwdim 2 -nldim 3 -sopt" + "./mixed_nonlinear_diffusion -online -tf 0.01 -rrdim 2 -rwdim 2 -nldim 3 -hrtype deim" + "./mixed_nonlinear_diffusion -online -tf 0.01 -rrdim 2 -rwdim 2 -nldim 3 -hrtype qdeim" + "./mixed_nonlinear_diffusion -online -tf 0.01 -rrdim 2 -rwdim 2 -nldim 3 -hrtype sopt" "./mixed_nonlinear_diffusion -online -tf 0.01 -rrdim 2 -rwdim 2 -nldim 3 -ns 1 -eqp -maxnnls 4" ) TYPE="PROM"