From c7edc3666ffdb0a7f4a15e42af04754e060f302d Mon Sep 17 00:00:00 2001 From: Chris Val Date: Tue, 27 Jun 2023 16:15:30 -0700 Subject: [PATCH 01/39] form combined NNLS problem --- rom/laghos_rom.cpp | 383 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 374 insertions(+), 9 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 7578f983..47cb6f0d 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -282,6 +282,7 @@ void ComputeElementRowOfG_V(const IntegrationRule *ir, trial_fe.CalcShape(ip, shape); // Compute the inner product, on this element, with the j-th V basis vector. + // TODO: Why compute shape * unitE just to get 1.0 back? r[q] = (v_e * Vloc_force) * (shape * unitE); } } @@ -338,6 +339,84 @@ void ComputeElementRowOfG_E(const IntegrationRule *ir, } } +// Sets the rows of constraints matrix G for the energy-conserving EQP rule. +void ComputeRowsOfG(const IntegrationRule *ir, + hydrodynamics::QuadratureData const& quad_data, + Vector const& w_e, Vector const& v_e, + FiniteElement const& test_fe, + FiniteElement const& trial_fe, + const int zone_id, + const bool equationV, const bool equationE, + Vector & rv, Vector & re) +{ + const int nqp = ir->GetNPoints(); + const int dim = trial_fe.GetDim(); // TODO: shouldn't it be the dimension of test_fe? + const int h1dofs_cnt = test_fe.GetDof(); + const int l2dofs_cnt = trial_fe.GetDof(); + + if (equationV) + { + MFEM_VERIFY(rv.Size() == nqp, ""); + } + + if (equationE) + { + MFEM_VERIFY(re.Size() == nqp, ""); + MFEM_VERIFY(v_e.Size() == h1dofs_cnt*dim, ""); + MFEM_VERIFY(w_e.Size() == l2dofs_cnt, ""); + } + + DenseMatrix grad_vshape(h1dofs_cnt, dim); // grad of velocity basis vector + DenseMatrix loc_force(h1dofs_cnt, dim); + Vector Vloc_force(loc_force.Data(), h1dofs_cnt * dim); + + for (int q = 0; q < nqp; q++) + { + const IntegrationPoint &ip = ir->IntPoint(q); + + // Form stress:grad_vshape at the current point. + // TODO: how do we set which velocity basis vector is used for the + // gradient calculation? + test_fe.CalcDShape(ip, grad_vshape); + for (int i = 0; i < h1dofs_cnt; i++) + { + for (int vd = 0; vd < dim; vd++) // Velocity components. + { + loc_force(i, vd) = 0.0; + for (int gd = 0; gd < dim; gd++) // Gradient components. + { + loc_force(i, vd) += + quad_data.stressJinvT(vd)(zone_id*nqp + q, gd) * grad_vshape(i,gd); + } + } + } + + // NOTE: UpdateQuadratureData includes ip.weight as a factor in quad_data.stressJinvT, + // set by LagrangianHydroOperator::UpdateQuadratureData. + loc_force *= 1.0 / ip.weight; // Divide by exact quadrature weight + + if (equationV) + { + // set rv = sum(Vloc_force) + rv[q] = 0.0; + for (i = 0; i < h1dofs_cnt * dim; i++) + { + rv[q] += Vloc_force(i); + } + } + if (equationE) + { + // set re = F^T * v; F = integral of (stress : grad w_i)phi_j + re[q] = 0.0; + for (i = 0; i < l2dofs_cnt; i++) + { + re[q] += w_e(i); + } + re[q] *= v_e * Vloc_force; + } + } // q -- quadrature point +} + #include "linalg/NNLS.h" void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, @@ -436,16 +515,34 @@ void ROM_Sampler::SetupEQP_Force_Eq(const CAROM::Matrix* snapX, Vector r(nqe); - // Compute G of size (NB * (nsnap+1)) x NQ, storing its transpose Gt. + // G is the matrix of accuracy constraints used to enforce that the + // evaluated quantity remains close to the result of the full quadrature + // case. + + // Compute G of size (NB * (nsnap+1)) x NQ, storing its transpose Gt. CAROM::Matrix Gt(NQ, NB * (nsnap+1), true); + // Velocity equation. // For 0 <= j < NB, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, - // G(j + (i*NB), (e*nqe) + m) - // is the coefficient of e_j^T M_e^{-1} F(v_i,e_i,x_i)^T v_i at point m of - // element e, with respect to the integration rule weight at that point, - // where the "exact" quadrature solution is ir0->GetWeights(). - - Vector v_i(tH1size); + // entry G(j + (i*NB), (e*nqe) + m) is the coefficient of + // + // e_j^T M_v^{-1} F(v_i,e_i,x_i)^T 1E + // + // at point m of element e with respect to the integration rule weight at + // that point, where the "exact" quadrature solution is ir0->GetWeights(). + // In the above, e_j is the jth velocity basis vector, 1E is the identity + // in the energy space. + + // Energy equation. + // Similarly, for 0 <= j < NB, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, + // entry G(j + (i*NB), (e*nqe) + m) is the coefficient of + // + // e_j^T M_e^{-1} F(v_i,e_i,x_i)^T v_i + // + // where e_j is the jth energy basis vector, v_i is the ith velocity + // snapshot. + + Vector v_i(tH1size); Vector x_i(tH1size); Vector e_i(tL2size); @@ -597,6 +694,267 @@ void ROM_Sampler::SetupEQP_Force_Eq(const CAROM::Matrix* snapX, std::to_string(rank)); } +void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, + const CAROM::Matrix* snapV, + const CAROM::Matrix* snapE, + const CAROM::Matrix* basisV, + const CAROM::Matrix* basisE, + ROM_Options const& input) +{ + const IntegrationRule *ir0 = input.FOMoper->GetIntegrationRule(); + const int nqe = ir0->GetNPoints(); + const int ne = input.H1FESpace->GetNE(); + const int NQ = ne * nqe; + + const int NBv = basisV->numColumns(); + const int NBe = basisE->numColumns(); + const int NBmin = min(NBv, NBe); + + Array numSnapVar(3); + numSnapVar[0] = snapX->numColumns(); + numSnapVar[1] = snapV->numColumns(); + numSnapVar[2] = snapE->numColumns(); + + const int nsnap = numSnapVar.Max(); + + Array numSkipped(3); + for (int i=0; i<3; ++i) numSkipped[i] = nsnap - numSnapVar[i]; + MFEM_VERIFY(numSkipped.Max() <= 1, ""); + + Vector rv(nqe), re(nqe); + + // G is the matrix of accuracy constraints used to enforce that the + // evaluated quantity remains close to the result of the full quadrature rule. + + // Declare G of size ((NBv + NBe) * (nsnap+1)) x NQ; store its transpose Gt. + // The first NBv*(nsnap+1) rows of G hold the velocity constraints; + // the remaining NBe*(nsnap+1) rows hold the energy constraints. + CAROM::Matrix Gt(NQ, (NBv + NBe) * (nsnap+1), true); + + // row index of G where energy constraints start + int estart = NBv * (nsnap + 1); + + Vector v_i(tH1size), x_i(tH1size), e_i(tL2size); + Vector w_j_e, v_i_e, v_j_e; + + Vector S((2*input.H1FESpace->GetVSize()) + input.L2FESpace->GetVSize()); + Vector S_v(S, input.H1FESpace->GetVSize(), input.H1FESpace->GetVSize()); // Subvector + + MFEM_VERIFY(tH1size == basisV->numRows(), ""); + MFEM_VERIFY(tL2size == basisE->numRows(), ""); + CAROM::Matrix Wv(H1size, NBv, true); + CAROM::Matrix We(L2size, NBe, true); + + // jth Wv columnn holds the jth velocity basis vector + for (int j=0; j const& w_el = ir0->GetWeights(); + MFEM_VERIFY(w_el.Size() == nqe, ""); + + for (int i=0; iResetQuadratureData(); + input.FOMoper->GetTimeStepEstimate(S); // Call this to call UpdateQuadratureData + input.FOMoper->ResetQuadratureData(); + + // Velocity equation. + // For 0 <= j < NBv, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, + // entry G(j + (i*NBv), (e*nqe) + m) is the coefficient of + // + // stress(v_i,e_i,x_i) : grad w_j + // + // at point m of element e with respect to the integration rule weight at + // that point, where the "exact" quadrature solution is ir0->GetWeights(). + // In the above, w_j is the jth velocity basis vector. + + // Energy equation. + // For 0 <= j < NBe, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, + // entry G(estart + j + (i*NBe), (e*nqe) + m) is the coefficient of + // + // (stress(v_i,e_i,x_i) : grad v_i) w_j + // + // at point m of element e with respect to the integration rule weight at + // that point, where the "exact" quadrature solution is ir0->GetWeights(). + // In the above, w_j is the jth energy basis vector. + + // Set the contraints for velocity and energy at the same time, then + // add the rest for velocity or energy, depending on which variable has + // more basis vectors (if not equal). + + for (int j=0; j Ediff(basisE->numRows()); + for (int k = 0; k < basisE->numRows(); ++k) + { + Ediff[k] = We(k, j); + } + gfL2.SetFromTrueDofs(Ediff); // jth energy basis vector + gfH1 = S_v; // ith velocity snapshot + + for (int e=0; eGetQuadData(), + w_j_e, v_i_e, + *input.H1FESpace->GetFE(e), + *input.L2FESpace->GetFE(e), + e, true, true, rv, re); + + for (int m=0; m NBmin = NBe, so there's more velocity constraints + for (int j=NBmin; jGetQuadData(), + w_j_e, v_i_e, + *input.H1FESpace->GetFE(e), + *input.L2FESpace->GetFE(e), + e, true, false, rv, re); + + for (int m=0; m NBmin = NBv, so there's more energy constraints + for (int j=NBmin; j Ediff(basisE->numRows()); + for (int k = 0; k < basisE->numRows(); ++k) + { + Ediff[k] = We(k, j); // jth energy basis vector + } + gfL2.SetFromTrueDofs(Ediff); + gfH1 = S_v; // ith velocity snapshot + + for (int e=0; eGetQuadData(), + w_j_e, v_i_e, + *input.H1FESpace->GetFE(e), + *input.L2FESpace->GetFE(e), + e, false, true, rv, re); + + for (int m=0; mnumRows() == input.H1FESpace->GetTrueVSize(), ""); MFEM_VERIFY(snapE->numRows() == input.L2FESpace->GetTrueVSize(), ""); - SetupEQP_Force_Eq(snapX, snapV, snapE, basisV, basisE, input, false); - SetupEQP_Force_Eq(snapX, snapV, snapE, basisV, basisE, input, true); + if (input.hyperreductionSamplingType == "eqp") + { + SetupEQP_Force_Eq(snapX, snapV, snapE, basisV, basisE, input, false); + SetupEQP_Force_Eq(snapX, snapV, snapE, basisV, basisE, input, true); + } + else if (input.hyperreductionSamplingType == "eqp_energy") + { + SetupEQP_En_Force_Eq(snapX, snapV, snapE, basisV, basisE, input); + } } void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) From b78829d3af41743a61bfaf74c68b21ee5df50640 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Wed, 28 Jun 2023 09:59:00 -0700 Subject: [PATCH 02/39] modify read/write of NNLS solution --- rom/laghos_rom.cpp | 430 ++++++++++++++++++++++++--------------------- rom/laghos_rom.hpp | 3 + 2 files changed, 228 insertions(+), 205 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 47cb6f0d..e6c4b70d 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -341,65 +341,65 @@ void ComputeElementRowOfG_E(const IntegrationRule *ir, // Sets the rows of constraints matrix G for the energy-conserving EQP rule. void ComputeRowsOfG(const IntegrationRule *ir, - hydrodynamics::QuadratureData const& quad_data, - Vector const& w_e, Vector const& v_e, - FiniteElement const& test_fe, - FiniteElement const& trial_fe, - const int zone_id, - const bool equationV, const bool equationE, - Vector & rv, Vector & re) + hydrodynamics::QuadratureData const& quad_data, + Vector const& w_e, Vector const& v_e, + FiniteElement const& test_fe, + FiniteElement const& trial_fe, + const int zone_id, + const bool equationV, const bool equationE, + Vector & rv, Vector & re) { - const int nqp = ir->GetNPoints(); - const int dim = trial_fe.GetDim(); // TODO: shouldn't it be the dimension of test_fe? - const int h1dofs_cnt = test_fe.GetDof(); - const int l2dofs_cnt = trial_fe.GetDof(); + const int nqp = ir->GetNPoints(); + const int dim = trial_fe.GetDim(); // TODO: shouldn't it be the dimension of test_fe? + const int h1dofs_cnt = test_fe.GetDof(); + const int l2dofs_cnt = trial_fe.GetDof(); if (equationV) { MFEM_VERIFY(rv.Size() == nqp, ""); } - + if (equationE) { MFEM_VERIFY(re.Size() == nqp, ""); MFEM_VERIFY(v_e.Size() == h1dofs_cnt*dim, ""); - MFEM_VERIFY(w_e.Size() == l2dofs_cnt, ""); + MFEM_VERIFY(w_e.Size() == l2dofs_cnt, ""); } - DenseMatrix grad_vshape(h1dofs_cnt, dim); // grad of velocity basis vector + DenseMatrix grad_vshape(h1dofs_cnt, dim); // grad of velocity basis vector DenseMatrix loc_force(h1dofs_cnt, dim); - Vector Vloc_force(loc_force.Data(), h1dofs_cnt * dim); + Vector Vloc_force(loc_force.Data(), h1dofs_cnt * dim); - for (int q = 0; q < nqp; q++) - { - const IntegrationPoint &ip = ir->IntPoint(q); + for (int q = 0; q < nqp; q++) + { + const IntegrationPoint &ip = ir->IntPoint(q); - // Form stress:grad_vshape at the current point. + // Form stress:grad_vshape at the current point. // TODO: how do we set which velocity basis vector is used for the // gradient calculation? - test_fe.CalcDShape(ip, grad_vshape); - for (int i = 0; i < h1dofs_cnt; i++) - { - for (int vd = 0; vd < dim; vd++) // Velocity components. - { - loc_force(i, vd) = 0.0; - for (int gd = 0; gd < dim; gd++) // Gradient components. - { - loc_force(i, vd) += - quad_data.stressJinvT(vd)(zone_id*nqp + q, gd) * grad_vshape(i,gd); - } - } - } + test_fe.CalcDShape(ip, grad_vshape); + for (int i = 0; i < h1dofs_cnt; i++) + { + for (int vd = 0; vd < dim; vd++) // Velocity components. + { + loc_force(i, vd) = 0.0; + for (int gd = 0; gd < dim; gd++) // Gradient components. + { + loc_force(i, vd) += + quad_data.stressJinvT(vd)(zone_id*nqp + q, gd) * grad_vshape(i,gd); + } + } + } - // NOTE: UpdateQuadratureData includes ip.weight as a factor in quad_data.stressJinvT, - // set by LagrangianHydroOperator::UpdateQuadratureData. - loc_force *= 1.0 / ip.weight; // Divide by exact quadrature weight + // NOTE: UpdateQuadratureData includes ip.weight as a factor in quad_data.stressJinvT, + // set by LagrangianHydroOperator::UpdateQuadratureData. + loc_force *= 1.0 / ip.weight; // Divide by exact quadrature weight if (equationV) { // set rv = sum(Vloc_force) rv[q] = 0.0; - for (i = 0; i < h1dofs_cnt * dim; i++) + for (int i = 0; i < h1dofs_cnt * dim; i++) { rv[q] += Vloc_force(i); } @@ -408,13 +408,13 @@ void ComputeRowsOfG(const IntegrationRule *ir, { // set re = F^T * v; F = integral of (stress : grad w_i)phi_j re[q] = 0.0; - for (i = 0; i < l2dofs_cnt; i++) + for (int i = 0; i < l2dofs_cnt; i++) { re[q] += w_e(i); } - re[q] *= v_e * Vloc_force; + re[q] *= v_e * Vloc_force; } - } // q -- quadrature point + } // q -- quadrature point } #include "linalg/NNLS.h" @@ -488,6 +488,7 @@ void WriteSolutionNNLS(CAROM::Vector const& sol, const string filename) outfile.close(); } +// Compute the reduced quadrature rules for the basic EQP case. void ROM_Sampler::SetupEQP_Force_Eq(const CAROM::Matrix* snapX, const CAROM::Matrix* snapV, const CAROM::Matrix* snapE, @@ -523,10 +524,10 @@ void ROM_Sampler::SetupEQP_Force_Eq(const CAROM::Matrix* snapX, CAROM::Matrix Gt(NQ, NB * (nsnap+1), true); // Velocity equation. - // For 0 <= j < NB, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, - // entry G(j + (i*NB), (e*nqe) + m) is the coefficient of + // For 0 <= j < NB, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, + // entry G(j + (i*NB), (e*nqe) + m) is the coefficient of // - // e_j^T M_v^{-1} F(v_i,e_i,x_i)^T 1E + // e_j^T M_v^{-1} F(v_i,e_i,x_i)^T 1E // // at point m of element e with respect to the integration rule weight at // that point, where the "exact" quadrature solution is ir0->GetWeights(). @@ -534,10 +535,10 @@ void ROM_Sampler::SetupEQP_Force_Eq(const CAROM::Matrix* snapX, // in the energy space. // Energy equation. - // Similarly, for 0 <= j < NB, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, - // entry G(j + (i*NB), (e*nqe) + m) is the coefficient of + // Similarly, for 0 <= j < NB, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, + // entry G(j + (i*NB), (e*nqe) + m) is the coefficient of // - // e_j^T M_e^{-1} F(v_i,e_i,x_i)^T v_i + // e_j^T M_e^{-1} F(v_i,e_i,x_i)^T v_i // // where e_j is the jth energy basis vector, v_i is the ith velocity // snapshot. @@ -694,58 +695,60 @@ void ROM_Sampler::SetupEQP_Force_Eq(const CAROM::Matrix* snapX, std::to_string(rank)); } +// Compute the reduced quadrature rule for the energy-conserving EQP case. void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, - const CAROM::Matrix* snapV, - const CAROM::Matrix* snapE, - const CAROM::Matrix* basisV, - const CAROM::Matrix* basisE, - ROM_Options const& input) + const CAROM::Matrix* snapV, + const CAROM::Matrix* snapE, + const CAROM::Matrix* basisV, + const CAROM::Matrix* basisE, + ROM_Options const& input) { - const IntegrationRule *ir0 = input.FOMoper->GetIntegrationRule(); - const int nqe = ir0->GetNPoints(); - const int ne = input.H1FESpace->GetNE(); - const int NQ = ne * nqe; - + const IntegrationRule *ir0 = input.FOMoper->GetIntegrationRule(); + const int nqe = ir0->GetNPoints(); + const int ne = input.H1FESpace->GetNE(); + const int NQ = ne * nqe; + const int NBv = basisV->numColumns(); const int NBe = basisE->numColumns(); const int NBmin = min(NBv, NBe); - Array numSnapVar(3); - numSnapVar[0] = snapX->numColumns(); - numSnapVar[1] = snapV->numColumns(); - numSnapVar[2] = snapE->numColumns(); + Array numSnapVar(3); + numSnapVar[0] = snapX->numColumns(); + numSnapVar[1] = snapV->numColumns(); + numSnapVar[2] = snapE->numColumns(); - const int nsnap = numSnapVar.Max(); + const int nsnap = numSnapVar.Max(); - Array numSkipped(3); - for (int i=0; i<3; ++i) numSkipped[i] = nsnap - numSnapVar[i]; - MFEM_VERIFY(numSkipped.Max() <= 1, ""); + Array numSkipped(3); + for (int i=0; i<3; ++i) numSkipped[i] = nsnap - numSnapVar[i]; + MFEM_VERIFY(numSkipped.Max() <= 1, ""); - Vector rv(nqe), re(nqe); + Vector rv(nqe), re(nqe); // G is the matrix of accuracy constraints used to enforce that the // evaluated quantity remains close to the result of the full quadrature rule. - + // Declare G of size ((NBv + NBe) * (nsnap+1)) x NQ; store its transpose Gt. // The first NBv*(nsnap+1) rows of G hold the velocity constraints; // the remaining NBe*(nsnap+1) rows hold the energy constraints. - CAROM::Matrix Gt(NQ, (NBv + NBe) * (nsnap+1), true); + CAROM::Matrix Gt(NQ, (NBv + NBe) * (nsnap+1), true); // row index of G where energy constraints start - int estart = NBv * (nsnap + 1); + const int estart = NBv * (nsnap + 1); - Vector v_i(tH1size), x_i(tH1size), e_i(tL2size); - Vector w_j_e, v_i_e, v_j_e; + Vector v_i(tH1size), x_i(tH1size), e_i(tL2size); + Vector w_j_e, v_i_e, v_j_e; - Vector S((2*input.H1FESpace->GetVSize()) + input.L2FESpace->GetVSize()); - Vector S_v(S, input.H1FESpace->GetVSize(), input.H1FESpace->GetVSize()); // Subvector + Vector S((2*input.H1FESpace->GetVSize()) + input.L2FESpace->GetVSize()); + Vector S_v(S, input.H1FESpace->GetVSize(), input.H1FESpace->GetVSize()); // Subvector - MFEM_VERIFY(tH1size == basisV->numRows(), ""); - MFEM_VERIFY(tL2size == basisE->numRows(), ""); + MFEM_VERIFY(tH1size == basisV->numRows(), ""); + MFEM_VERIFY(tL2size == basisE->numRows(), ""); CAROM::Matrix Wv(H1size, NBv, true); CAROM::Matrix We(L2size, NBe, true); // jth Wv columnn holds the jth velocity basis vector + // TODO: do we need Wv anymore? it's not used anywhere for (int j=0; j const& w_el = ir0->GetWeights(); - MFEM_VERIFY(w_el.Size() == nqe, ""); + Array const& w_el = ir0->GetWeights(); + MFEM_VERIFY(w_el.Size() == nqe, ""); - for (int i=0; iResetQuadratureData(); - input.FOMoper->GetTimeStepEstimate(S); // Call this to call UpdateQuadratureData - input.FOMoper->ResetQuadratureData(); + // NOTE: after SetStateFromTrueDOFs, gfH1 is the V-component of S + input.FOMoper->ResetQuadratureData(); + input.FOMoper->GetTimeStepEstimate(S); // to call UpdateQuadratureData + input.FOMoper->ResetQuadratureData(); // Velocity equation. // For 0 <= j < NBv, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, - // entry G(j + (i*NBv), (e*nqe) + m) is the coefficient of + // entry G(j + (i*NBv), (e*nqe) + m) is the coefficient of // - // stress(v_i,e_i,x_i) : grad w_j + // stress(v_i,e_i,x_i) : grad w_j // // at point m of element e with respect to the integration rule weight at // that point, where the "exact" quadrature solution is ir0->GetWeights(). // In the above, w_j is the jth velocity basis vector. // Energy equation. - // For 0 <= j < NBe, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, - // entry G(estart + j + (i*NBe), (e*nqe) + m) is the coefficient of + // For 0 <= j < NBe, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, + // entry G(estart + j + (i*NBe), (e*nqe) + m) is the coefficient of // // (stress(v_i,e_i,x_i) : grad v_i) w_j // // at point m of element e with respect to the integration rule weight at // that point, where the "exact" quadrature solution is ir0->GetWeights(). // In the above, w_j is the jth energy basis vector. - + // Set the contraints for velocity and energy at the same time, then // add the rest for velocity or energy, depending on which variable has // more basis vectors (if not equal). for (int j=0; jnumRows()); + for (int k = 0; k < basisE->numRows(); ++k) + { + Ediff[k] = We(k, j); + } + gfL2.SetFromTrueDofs(Ediff); // jth energy basis vector + gfH1 = S_v; // ith velocity snapshot - Array Ediff(basisE->numRows()); - for (int k = 0; k < basisE->numRows(); ++k) - { - Ediff[k] = We(k, j); - } - gfL2.SetFromTrueDofs(Ediff); // jth energy basis vector - gfH1 = S_v; // ith velocity snapshot - for (int e=0; eGetQuadData(), - w_j_e, v_i_e, - *input.H1FESpace->GetFE(e), - *input.L2FESpace->GetFE(e), - e, true, true, rv, re); + w_j_e, v_i_e, + *input.H1FESpace->GetFE(e), + *input.L2FESpace->GetFE(e), + e, true, true, rv, re); - for (int m=0; m NBmin = NBe, so there's more velocity constraints for (int j=NBmin; jGetQuadData(), - w_j_e, v_i_e, - *input.H1FESpace->GetFE(e), - *input.L2FESpace->GetFE(e), - e, true, false, rv, re); - - for (int m=0; mGetQuadData(), + w_j_e, v_i_e, + *input.H1FESpace->GetFE(e), + *input.L2FESpace->GetFE(e), + e, true, false, rv, re); + + for (int m=0; m NBmin = NBv, so there's more energy constraints for (int j=NBmin; j Ediff(basisE->numRows()); + Vector Ediff(basisE->numRows()); for (int k = 0; k < basisE->numRows(); ++k) { Ediff[k] = We(k, j); // jth energy basis vector } gfL2.SetFromTrueDofs(Ediff); gfH1 = S_v; // ith velocity snapshot - + for (int e=0; eGetQuadData(), - w_j_e, v_i_e, - *input.H1FESpace->GetFE(e), - *input.L2FESpace->GetFE(e), - e, false, true, rv, re); + w_j_e, v_i_e, + *input.H1FESpace->GetFE(e), + *input.L2FESpace->GetFE(e), + e, false, true, rv, re); for (int m=0; mnumRows() == input.H1FESpace->GetTrueVSize(), ""); MFEM_VERIFY(basisE->numRows() == input.L2FESpace->GetTrueVSize(), ""); @@ -966,13 +971,15 @@ void ROM_Sampler::SetupEQP_Force(const CAROM::Matrix* snapX, const CAROM::Matrix MFEM_VERIFY(snapV->numRows() == input.H1FESpace->GetTrueVSize(), ""); MFEM_VERIFY(snapE->numRows() == input.L2FESpace->GetTrueVSize(), ""); - if (input.hyperreductionSamplingType == "eqp") + if (input.hyperreductionSamplingType == eqp) { + // basic EQP: different rules for velocity and energy SetupEQP_Force_Eq(snapX, snapV, snapE, basisV, basisE, input, false); SetupEQP_Force_Eq(snapX, snapV, snapE, basisV, basisE, input, true); } - else if (input.hyperreductionSamplingType == "eqp_energy") + else if (input.hyperreductionSamplingType == eqp_energy) { + // energy-conserving EQP: one combined rule for velocity and energy SetupEQP_En_Force_Eq(snapX, snapV, snapE, basisV, basisE, input); } } @@ -2724,16 +2731,29 @@ ROM_Operator::ROM_Operator(ROM_Options const& input, ROM_Basis *b, ComputeReducedMe(); } - if (hyperreduce && hyperreductionSamplingType == eqp) - { - // TODO: are reduced mass matrices needed for EQP, or just W matrices? + if (hyperreduce && hyperreductionSamplingType == eqp) + { + // basic EQP + // TODO: are reduced mass matrices needed for EQP, or just W matrices? - ComputeReducedMv(); - ComputeReducedMe(); + // TODO: why compute the reduced mass matrices, when for the basic EQP + // case we've already multiplied my M^{-1} and basis^T when deriving + // the reduced quadrature rules? - ReadSolutionNNLS(input, "run/nnlsV", eqpI, eqpW); - ReadSolutionNNLS(input, "run/nnlsE", eqpI_E, eqpW_E); - } + // ComputeReducedMv(); + // ComputeReducedMe(); + + ReadSolutionNNLS(input, "run/nnlsV", eqpI, eqpW); + ReadSolutionNNLS(input, "run/nnlsE", eqpI_E, eqpW_E); + } + else if (hyperreduce && hyperreductionSamplingType == eqp_energy) + { + // energy-conserving EQP + // NOTE: will not compute reduced M matrices; instead, will normalize + // the basis matrices so that the reduced M matrices are identity. + + ReadSolutionNNLS(input, "run/nnlsEC", eqpI, eqpW); + } } void ROM_Operator::ReadSolutionNNLS(ROM_Options const& input, string basename, diff --git a/rom/laghos_rom.hpp b/rom/laghos_rom.hpp index 3104bd50..3619d2f5 100644 --- a/rom/laghos_rom.hpp +++ b/rom/laghos_rom.hpp @@ -822,6 +822,9 @@ class ROM_Sampler void SetupEQP_Force_Eq(const CAROM::Matrix* snapX, const CAROM::Matrix* snapV, const CAROM::Matrix* snapE, const CAROM::Matrix* basisV, const CAROM::Matrix* basisE, ROM_Options const& input, bool equationE); + + void SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, const CAROM::Matrix* snapV, const CAROM::Matrix* snapE, + const CAROM::Matrix* basisV, const CAROM::Matrix* basisE, ROM_Options const& input); }; class ROM_Basis From c4280afc19aaf0f1a74de2014466b6311a0c1c13 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Fri, 30 Jun 2023 14:11:00 -0700 Subject: [PATCH 03/39] form RHS force vectors --- rom/laghos.cpp | 10 ++- rom/laghos_rom.cpp | 167 +++++++++++++++++++++++++----------------- rom/laghos_rom.hpp | 16 ++-- rom/laghos_solver.cpp | 8 +- 4 files changed, 124 insertions(+), 77 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 1aa8c8f9..a0e36ea8 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -1573,7 +1573,15 @@ int main(int argc, char *argv[]) { romOper[0]->ApplyHyperreduction(romS); } - double windowEndpoint = 0.0; + + // To perform the induced Gram-Schmidt, so that the reduced mass + // matrices are identity. + if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) + { + romOper[0]->ApplyHyperreduction(romS); + } + + double windowEndpoint = 0.0; double windowOverlapMidpoint = 0.0; for (ti = 1; !last_step; ti++) { diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index e6c4b70d..13170743 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -2084,6 +2084,9 @@ void ROM_Basis::ComputeReducedMatrices(bool sns1) MFEM_VERIFY(BX0->dim() == rdimx, ""); } + // TODO: what do this and the following if-blocks compute? + // Which of those computations are needed in the energy-conserving + // EQP? if (!sns1) { // Compute reduced matrix BsinvV = (BVsp^T BFvsp BsinvV^T)^T = BsinvV BFvsp^T BVsp @@ -2105,6 +2108,8 @@ void ROM_Basis::ComputeReducedMatrices(bool sns1) BsinvE = prod2; } + // TODO: We are using RK2-avg in energy-conserving EQP; do we need + // the matrices computed here though? if (RK2AvgFormulation) { const CAROM::Matrix *prodX = BXsp->transposeMult(BXsp); @@ -3412,6 +3417,7 @@ void ROM_Operator::UndoInducedGramSchmidt(const int var, Vector &S, bool keep_da void ROM_Operator::ApplyHyperreduction(Vector &S) { + // TODO: should set GramSchmidt = true when using energy-conserving EQP. if (useGramSchmidt && !sns1) { InducedGramSchmidt(1, S); // velocity @@ -4397,7 +4403,8 @@ CAROM::GreedySampler* UseROMDatabase(ROM_Options& romOptions, const int myid, co return parameterPointGreedySampler; } -void ROM_Operator::ForceIntegratorEQP(Vector & res) const +void ROM_Operator::ForceIntegratorEQP(Vector & res, + bool energy_conserve) const { const IntegrationRule *ir = operFOM->GetIntegrationRule(); const int rdim = basis->GetDimV(); @@ -4407,13 +4414,9 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res) const const int nqe = ir->GetWeights().Size(); DenseMatrix vshape, loc_force; - Vector shape, unitE, rhs; - Array vdofs; - - Vector v_e; - - res = 0.0; + + res = 0.0; int eprev = -1; int dof = 0; @@ -4501,15 +4504,13 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res) const // TODO: reduce this UpdateQuadratureData function to the EQP points. const int h1dofs_cnt = test_fe->GetDof(); - const int dim = trial_fe->GetDim(); + const int dim = trial_fe->GetDim(); // TODO: shouldn't it be the dim of test_fe? const int l2dofs_cnt = trial_fe->GetDof(); vshape.SetSize(h1dofs_cnt, dim); loc_force.SetSize(h1dofs_cnt, dim); - shape.SetSize(l2dofs_cnt); - // Form stress:grad_shape at the current point. test_fe->CalcDShape(ip, vshape); for (int k = 0; k < h1dofs_cnt; k++) @@ -4524,36 +4525,52 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res) const } } } - loc_force *= eqpW[i] / ip.weight; // Replace exact quadrature weight with EQP weight. + + Vector Vloc_force(loc_force.Data(), loc_force.NumRows() * loc_force.NumCols()); + Vector v_e(h1dofs_cnt * dim); - trial_fe->CalcShape(ip, shape); - - Vector Vloc_force(loc_force.Data(), loc_force.NumRows() * loc_force.NumCols()); - - unitE.SetSize(shape.Size()); - unitE = 1.0; + if (energy_conserve) + { + // energy-conserving EQP + for (int j = 0; j < rdim; ++j) + { + // TODO: v_e should be the jth velocity basis vector DOFs that + // correspond to the current element. + // Recall that the basis vectors have been renormalized by the + // induced Gram-Schmidt process, so that reduced Mv = I. - rhs.SetSize(h1dofs_cnt * dim); + // I think the below is not right. + basis->GetBasisVectorV(false, j, v_e); - v_e.SetSize(rhs.Size()); + res[j] += v_e * Vloc_force; + } + } + else + { + // basic EQP + Vector shape(l2dofs_cnt), unitE(l2dofs_cnt); + trial_fe->CalcShape(ip, shape); + unitE = 1.0; - const int eos = elemIndex * nvdof; - for (int j=0; jGetIntegrationRule(); const int rdim = basis->GetDimE(); @@ -4563,11 +4580,9 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res) const const int nqe = ir->GetWeights().Size(); DenseMatrix vshape, loc_force; - Vector shape, rhs; - Array vdofs, edofs; - Vector v_e, w_e; + Vector v_e; res = 0.0; @@ -4679,8 +4694,6 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res) const vshape.SetSize(h1dofs_cnt, dim); loc_force.SetSize(h1dofs_cnt, dim); - shape.SetSize(l2dofs_cnt); - // Form stress:grad_shape at the current point. test_fe->CalcDShape(ip, vshape); for (int k = 0; k < h1dofs_cnt; k++) @@ -4696,40 +4709,54 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res) const } } - loc_force *= eqpW_E[i] / ip.weight; // Replace exact quadrature weight with EQP weight. + loc_force *= eqpW_E[i] / ip.weight; // Replace exact quadrature weight with EQP weight. - trial_fe->CalcShape(ip, shape); + Vector Vloc_force(loc_force.Data(), loc_force.NumRows() * loc_force.NumCols()); - Vector Vloc_force(loc_force.Data(), loc_force.NumRows() * loc_force.NumCols()); + Vector shape(l2dofs_cnt); + trial_fe->CalcShape(ip, shape); - rhs.SetSize(h1dofs_cnt * dim); + if (energy_conserve) + { + // energy-conserving EQP + for (int j = 0; j < rdim; j++) + { + for (int k = 0; k < l2dofs_cnt; k++) res[j] += shape(k); - w_e.SetSize(nedof); + res[j] *= v_e * Vloc_force; + } + } + else + { + // basic EQP + Vector w_e(nedof); + const int eos = elemIndex * nedof; - const int eos = elemIndex * nedof; - for (int j=0; jGetDimV()); - ForceIntegratorEQP(res); + ForceIntegratorEQP(res, energy_conserve); basis->LiftROMtoFOM_dVdt(res, rhs); } -void ROM_Operator::ForceIntegratorEQP_E_FOM(Vector const& v, Vector & rhs) const +void ROM_Operator::ForceIntegratorEQP_E_FOM(Vector const& v, Vector & rhs, + bool energy_conserve) const { Vector res(basis->GetDimE()); @@ -4750,10 +4777,9 @@ void ROM_Operator::StepRK2Avg(Vector &S, double &t, double &dt) const else basis->LiftROMtoFOM(S, fx); - if (hyperreduce && hyperreductionSamplingType == eqp) - { - operFOM->SetRomOperator(this); - } + if (hyperreduce) + if (hyperreductionSamplingType == eqp || hyperreductionSamplingType == eqp_energy) + operFOM->SetRomOperator(this); const int Vsize = use_sample_mesh ? basis->SolutionSizeH1SP() : basis->SolutionSizeH1FOM(); const int Esize = basis->SolutionSizeL2SP(); @@ -4773,36 +4799,41 @@ void ROM_Operator::StepRK2Avg(Vector &S, double &t, double &dt) const // - Update V using dv_dt. // - Compute de_dt and dx_dt using S and V. - // -- 1. - // S is S0. + // -- Stage 1 of 2 (S is S_n = S0). hydro_oper->UpdateMesh(fx); + hydro_oper->SolveVelocity(fx, dS_dt); - if (use_sample_mesh) basis->HyperreduceRHS_V(dv_dt); // Set dv_dt based on RHS computed by SolveVelocity + if (use_sample_mesh) basis->HyperreduceRHS_V(dv_dt); // Set dv_dt based on RHS computed by SolveVelocity - // V = v0 + 0.5 * dt * dv_dt; + // time march velocity vector: V = V_{n+1/2} = v0 + 0.5 * dt * dv_dt add(v0, 0.5 * dt, dv_dt, V); + hydro_oper->SolveEnergy(fx, V, dS_dt); if (use_sample_mesh) basis->HyperreduceRHS_E(de_dt); // Set de_dt based on RHS computed by SolveEnergy - dx_dt = V; + dx_dt = V; - // -- 2. - // S = S0 + 0.5 * dt * dS_dt; + // time march full state vector: S = S_{n+1/2} = S0 + 0.5 * dt * dS_dt add(S0, 0.5 * dt, dS_dt, fx); + + // -- Stage 2 of 2 (S is S_{n+1/2}). hydro_oper->ResetQuadratureData(); hydro_oper->UpdateMesh(fx); + hydro_oper->SolveVelocity(fx, dS_dt); if (use_sample_mesh) basis->HyperreduceRHS_V(dv_dt); // Set dv_dt based on RHS computed by SolveVelocity - // V = v0 + 0.5 * dt * dv_dt; - add(v0, 0.5 * dt, dv_dt, V); - hydro_oper->SolveEnergy(fx, V, dS_dt); + + // time march velocity vector: V = v0 + 0.5 * dt * dv_dt + add(v0, 0.5 * dt, dv_dt, V); + + // V = average V_{n+1/2} + hydro_oper->SolveEnergy(fx, V, dS_dt); if (use_sample_mesh) basis->HyperreduceRHS_E(de_dt); // Set de_dt based on RHS computed by SolveEnergy - dx_dt = V; + dx_dt = V; - // -- 3. - // S = S0 + dt * dS_dt. + // time march full state vector: S = S_{n+1} = S0 + dt * dS_dt add(S0, dt, dS_dt, fx); - hydro_oper->ResetQuadratureData(); - + + hydro_oper->ResetQuadratureData(); MFEM_VERIFY(!useReducedM, "TODO"); if (use_sample_mesh) diff --git a/rom/laghos_rom.hpp b/rom/laghos_rom.hpp index 3619d2f5..8cf67e88 100644 --- a/rom/laghos_rom.hpp +++ b/rom/laghos_rom.hpp @@ -73,7 +73,8 @@ enum HyperreductionSamplingType gnat, // Default, GNAT qdeim, // QDEIM sopt, // S-OPT - eqp // EQP + eqp, // EQP + eqp_energy // Energy-conserving EQP }; static HyperreductionSamplingType getHyperreductionSamplingType(const char* sampling_type) @@ -83,7 +84,8 @@ static HyperreductionSamplingType getHyperreductionSamplingType(const char* samp {"gnat", gnat}, {"qdeim", qdeim}, {"sopt", sopt}, - {"eqp", eqp} + {"eqp", eqp}, + {"eqp_energy", eqp_energy} }; auto iter = SamplingTypeMap.find(sampling_type); MFEM_VERIFY(iter != std::end(SamplingTypeMap), "Invalid input for hyperreduction sampling type"); @@ -1225,11 +1227,13 @@ class ROM_Operator : public TimeDependentOperator void SolveSpaceTime(Vector &S); void SolveSpaceTimeGN(Vector &S); - void ForceIntegratorEQP_FOM(Vector & rhs) const; - void ForceIntegratorEQP(Vector & res) const; + void ForceIntegratorEQP_FOM(Vector & rhs, bool energy_conserve = false) const; + void ForceIntegratorEQP(Vector & res, bool energy_conserve = false) const; - void ForceIntegratorEQP_E_FOM(Vector const& v, Vector & rhs) const; - void ForceIntegratorEQP_E(Vector const& v, Vector & res) const; + void ForceIntegratorEQP_E_FOM(Vector const& v, Vector & rhs, + bool energy_conserve = false) const; + void ForceIntegratorEQP_E(Vector const& v, Vector & res, + bool energy_conserve = false) const; ~ROM_Operator() { diff --git a/rom/laghos_solver.cpp b/rom/laghos_solver.cpp index 2d433b2d..2773295b 100644 --- a/rom/laghos_solver.cpp +++ b/rom/laghos_solver.cpp @@ -277,7 +277,9 @@ void LagrangianHydroOperator::SolveVelocity(const Vector &S, if (use_eqp) { - rom_op->ForceIntegratorEQP_FOM(rhs); + // TODO: this function doesn't have access to the hyperreduction + // sampling type; how should I resolve this? + rom_op->ForceIntegratorEQP_FOM(rhs, hyperreduction.SamplingType == eqp_energy); } else { @@ -408,7 +410,9 @@ void LagrangianHydroOperator::SolveEnergy(const Vector &S, const Vector &v, timer.sw_force.Start(); if (use_eqp) { - rom_op->ForceIntegratorEQP_E_FOM(v, e_rhs); + // TODO: this function doesn't have access to the hyperreduction + // sampling type; how should I resolve this? + rom_op->ForceIntegratorEQP_E_FOM(v, e_rhs, hyperreduction.SamplingType == eqp_energy); } else { From 0da7e931ceaeb17ac451debb9497949b661c64c3 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Fri, 30 Jun 2023 14:50:39 -0700 Subject: [PATCH 04/39] modify parameter flags --- rom/laghos.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index a0e36ea8..fc399d37 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -394,8 +394,15 @@ int main(int argc, char *argv[]) romOptions.basisIdentifier = std::string(basisIdentifier); romOptions.hyperreductionSamplingType = getHyperreductionSamplingType(hyperreductionSamplingType); - romOptions.use_sample_mesh = romOptions.hyperreduce && (romOptions.hyperreductionSamplingType != eqp); - MFEM_VERIFY(!(romOptions.SNS) || (romOptions.hyperreductionSamplingType != eqp), "Using SNS with EQP is prohibited"); + romOptions.use_sample_mesh = romOptions.hyperreduce && (romOptions.hyperreductionSamplingType != eqp + || romOptions.hyperreductionSamplingType != eqp_energy); + + MFEM_VERIFY(!(romOptions.SNS) || (romOptions.hyperreductionSamplingType != eqp && + romOptions.hyperreductionSamplingType != eqp_energy), + "Using SNS with EQP is prohibited"); + + if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) + romOptions.GramSchmidt = true; romOptions.spaceTimeMethod = getSpaceTimeMethod(spaceTimeMethod); const bool spaceTime = (romOptions.spaceTimeMethod != no_space_time); From 90327aff65b9ea5db8f2b54e90eac39e5df4d594 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Thu, 6 Jul 2023 09:56:46 -0700 Subject: [PATCH 05/39] add energy identity to basis --- rom/laghos.cpp | 14 +++++++++--- rom/laghos_rom.cpp | 52 +++++++++++++++++++++++++++++++++++++++++++- rom/laghos_rom.hpp | 2 +- rom/laghos_utils.cpp | 8 +++++++ 4 files changed, 71 insertions(+), 5 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index fc399d37..b1e2d437 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -1581,9 +1581,9 @@ int main(int argc, char *argv[]) romOper[0]->ApplyHyperreduction(romS); } - // To perform the induced Gram-Schmidt, so that the reduced mass - // matrices are identity. - if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) + // To perform the induced Gram-Schmidt during the online stage, so that + // the reduced mass matrices are identity. + if (rom_online && romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) { romOper[0]->ApplyHyperreduction(romS); } @@ -1928,6 +1928,7 @@ int main(int argc, char *argv[]) } else { + // Form the reduced bases from the snapshots. sampler->Finalize(cutoff, romOptions); } if (myid == 0 && romOptions.parameterID == -1) { @@ -2057,6 +2058,13 @@ int main(int argc, char *argv[]) { romOper[romOptions.window]->ApplyHyperreduction(romS); } + + // To apply the induced Gram-Schmid to the basis of the + // new time window. + if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) + { + romOper[romOptions.window]->ApplyHyperreduction(romS); + } if (problem == 7 && romOptions.indicatorType == penetrationDistance) { diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 13170743..a72f7c53 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -862,6 +862,8 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, gfH1.GetElementDofValues(e, v_i_e); // set the constraints for the current basis vector & element + // rv: velocity constraints + // re: energy constraints ComputeRowsOfG(ir0, input.FOMoper->GetQuadData(), w_j_e, v_i_e, *input.H1FESpace->GetFE(e), @@ -986,6 +988,16 @@ void ROM_Sampler::SetupEQP_Force(const CAROM::Matrix* snapX, void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) { + if (input.hyperreductionSamplingType == eqp_energy) + { + // For the energy-conserving EQP case, we increase the energy basis + // dimension by 1 to accomodate for the addition of the energy + // identity. + // The change is done here, before the window basis parameters are + // written to their files, and is also carried to the caller's scope. + cutoff[2] += 1; + } + if (writeSnapshots) { if (!useXV) generator_X->writeSnapshot(); @@ -1091,6 +1103,44 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) delete tBasisE; } + if (input.hyperreductionSamplingType == eqp_energy) + { + const CAROM::Matrix *basisV = generator_V->getSpatialBasis(); + const CAROM::Matrix *basisE = generator_E->getSpatialBasis(); + + MPI_Bcast(cutoff.GetData(), cutoff.Size(), MPI_INT, 0, MPI_COMM_WORLD); + + // Truncate the bases. + // For the energy basis, the cutoff parameter has already been + // increased by 1 to accomodate the energy identity. + CAROM::Matrix *tBasisV = basisV->getFirstNColumns(cutoff[1]); + CAROM::Matrix *tBasisE = basisE->getFirstNColumns(cutoff[2]); + + // Form the energy identity and replace the last basis vector by it. + Vector unitE(input.tL2size); + + // TODO: my understanding is that the "=" opetaror is overloaded, so + // that, by setting unitE = 1.0, the vector unitE contains the coeffs + // that result in the summed function being 1.0 on the grid points. + // Is that correct? + unitE = 1.0; + + for (int i = 0; i < input.tL2size; i++) + { + // TODO: I believe I can access the last column in this way, + // without having to use pointer arithmetic, right? + (*tBasisE)(i, cutoff[2]-1) = unitE[i]; + } + + SetupEQP_Force(generator_X->getSnapshotMatrix(), + generator_V->getSnapshotMatrix(), + generator_E->getSnapshotMatrix(), + tBasisV, tBasisE, input); + + delete tBasisV; + delete tBasisE; + } + delete generator_X; delete generator_V; delete generator_E; @@ -2222,7 +2272,7 @@ void ROM_Basis::ReadSolutionBases(const int window) basisV = basisX; } - if (hyperreductionSamplingType == eqp) return; + if (hyperreductionSamplingType == eqp || hyperreductionSamplingType == eqp_energy) return; if (use_sns) // TODO: only do in online and not hyperreduce { diff --git a/rom/laghos_rom.hpp b/rom/laghos_rom.hpp index 8cf67e88..c79f46ba 100644 --- a/rom/laghos_rom.hpp +++ b/rom/laghos_rom.hpp @@ -481,7 +481,7 @@ class ROM_Sampler energyFraction_X(input.energyFraction_X), sns(input.SNS), lhoper(input.FOMoper), writeSnapshots(input.parameterID >= 0), parameterID(input.parameterID), basename(*input.basename), Voffset(!input.useXV && !input.useVX && !input.mergeXV), useXV(input.useXV), useVX(input.useVX), VTos(input.VTos), spaceTime(input.spaceTimeMethod != no_space_time), - rhsBasis(input.hyperreductionSamplingType != eqp) + rhsBasis(input.hyperreductionSamplingType != eqp && input.hyperreductionSamplingType != eqp_energy) { const int window = input.window; diff --git a/rom/laghos_utils.cpp b/rom/laghos_utils.cpp index fc5b14e4..8ac0e232 100644 --- a/rom/laghos_utils.cpp +++ b/rom/laghos_utils.cpp @@ -507,11 +507,19 @@ void ReadPDweight(std::vector& pd_weight, std::string outputPath) void SetWindowParameters(Array2D const& twparam, ROM_Options & romOptions) { const int w = romOptions.window; + romOptions.dimX = min(romOptions.max_dimX, twparam(w,0)); romOptions.dimV = min(romOptions.max_dimV, twparam(w,1)); + + // TODO: for the energy-conserving EQP case, if twparam(w,2) exceeds the + // set maximum, truncating the stored basis will lead to loss of the + // energy identity vector, thus violating one of the energy conservation + // conditions. romOptions.dimE = min(romOptions.max_dimE, twparam(w,2)); + romOptions.dimFv = romOptions.SNS ? romOptions.dimV : min(romOptions.max_dimFv, twparam(w,3)); romOptions.dimFe = romOptions.SNS ? romOptions.dimE : min(romOptions.max_dimFe, twparam(w,4)); + if (romOptions.useXV) romOptions.dimX = romOptions.dimV; if (romOptions.useVX) romOptions.dimV = romOptions.dimX; From 0b3907e3fc4f6789c00c5490b7c30b833ee3c2f7 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Fri, 7 Jul 2023 13:54:31 -0700 Subject: [PATCH 06/39] review 1 changes --- rom/laghos_rom.cpp | 9 +++++++-- rom/laghos_rom.hpp | 2 ++ rom/laghos_solver.cpp | 4 ++-- 3 files changed, 11 insertions(+), 4 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index a72f7c53..a54e7399 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -1117,7 +1117,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) CAROM::Matrix *tBasisE = basisE->getFirstNColumns(cutoff[2]); // Form the energy identity and replace the last basis vector by it. - Vector unitE(input.tL2size); + Vector unitE(tL2size); // TODO: my understanding is that the "=" opetaror is overloaded, so // that, by setting unitE = 1.0, the vector unitE contains the coeffs @@ -1125,7 +1125,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) // Is that correct? unitE = 1.0; - for (int i = 0; i < input.tL2size; i++) + for (int i = 0; i < tL2size; i++) { // TODO: I believe I can access the last column in this way, // without having to use pointer arithmetic, right? @@ -4899,3 +4899,8 @@ void ROM_Operator::StepRK2Avg(Vector &S, double &t, double &dt) const t += dt; } + +HyperreductionSamplingType ROM_Operator::getSamplingType() const +{ + return hyperreductionSamplingType; +} diff --git a/rom/laghos_rom.hpp b/rom/laghos_rom.hpp index c79f46ba..8a11719f 100644 --- a/rom/laghos_rom.hpp +++ b/rom/laghos_rom.hpp @@ -1235,6 +1235,8 @@ class ROM_Operator : public TimeDependentOperator void ForceIntegratorEQP_E(Vector const& v, Vector & res, bool energy_conserve = false) const; + HyperreductionSamplingType getSamplingType() const; + ~ROM_Operator() { delete mat_gf_coeff; diff --git a/rom/laghos_solver.cpp b/rom/laghos_solver.cpp index 2773295b..3b144f01 100644 --- a/rom/laghos_solver.cpp +++ b/rom/laghos_solver.cpp @@ -279,7 +279,7 @@ void LagrangianHydroOperator::SolveVelocity(const Vector &S, { // TODO: this function doesn't have access to the hyperreduction // sampling type; how should I resolve this? - rom_op->ForceIntegratorEQP_FOM(rhs, hyperreduction.SamplingType == eqp_energy); + rom_op->ForceIntegratorEQP_FOM(rhs, rom_op->getSamplingType() == eqp_energy); } else { @@ -412,7 +412,7 @@ void LagrangianHydroOperator::SolveEnergy(const Vector &S, const Vector &v, { // TODO: this function doesn't have access to the hyperreduction // sampling type; how should I resolve this? - rom_op->ForceIntegratorEQP_E_FOM(v, e_rhs, hyperreduction.SamplingType == eqp_energy); + rom_op->ForceIntegratorEQP_E_FOM(v, e_rhs, rom_op->getSamplingType() == eqp_energy); } else { From 9b9b0a75cdbac928d448fb284e75dfec934cfde3 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Tue, 11 Jul 2023 08:51:05 -0700 Subject: [PATCH 07/39] print energy total & diff --- rom/laghos.cpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index b1e2d437..9f1b2f92 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -2111,6 +2111,22 @@ int main(int argc, char *argv[]) << ",\t|e| = " << setprecision(10) << sqrt(tot_norm) << endl; } + + if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) + { + double energy_total, energy_diff; + energy_total = oper->InternalEnergy(*e_gf) + + oper->KineticEnergy(*v_gf); + energy_diff = energy_total - energy_init; + + if (mpi.Root()) + { + cout << "\tE_tot = " << setprecision(14) + << energy_total + << ",\tE_diff = " << setprecision(14) + << energy_diff << endl; + } + } // Make sure all ranks have sent their 'v' solution before initiating // another set of GLVis connections (one from each rank): From 3d661d2098a5cd32d9f358da63423214e6a26042 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Tue, 11 Jul 2023 15:11:42 -0700 Subject: [PATCH 08/39] review 2 changes --- rom/laghos.cpp | 7 +- rom/laghos_rom.cpp | 165 ++++++++++++++++++++++-------------------- rom/laghos_solver.cpp | 4 - 3 files changed, 88 insertions(+), 88 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 9f1b2f92..2a100947 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -402,7 +402,7 @@ int main(int argc, char *argv[]) "Using SNS with EQP is prohibited"); if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) - romOptions.GramSchmidt = true; + romOptions.GramSchmidt = false; romOptions.spaceTimeMethod = getSpaceTimeMethod(spaceTimeMethod); const bool spaceTime = (romOptions.spaceTimeMethod != no_space_time); @@ -1581,8 +1581,7 @@ int main(int argc, char *argv[]) romOper[0]->ApplyHyperreduction(romS); } - // To perform the induced Gram-Schmidt during the online stage, so that - // the reduced mass matrices are identity. + // TODO: do we want that for the energy-conserving EQP? if (rom_online && romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) { romOper[0]->ApplyHyperreduction(romS); @@ -2059,8 +2058,6 @@ int main(int argc, char *argv[]) romOper[romOptions.window]->ApplyHyperreduction(romS); } - // To apply the induced Gram-Schmid to the basis of the - // new time window. if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) { romOper[romOptions.window]->ApplyHyperreduction(romS); diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index a54e7399..6453a56e 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -282,7 +282,6 @@ void ComputeElementRowOfG_V(const IntegrationRule *ir, trial_fe.CalcShape(ip, shape); // Compute the inner product, on this element, with the j-th V basis vector. - // TODO: Why compute shape * unitE just to get 1.0 back? r[q] = (v_e * Vloc_force) * (shape * unitE); } } @@ -339,10 +338,11 @@ void ComputeElementRowOfG_E(const IntegrationRule *ir, } } +// TODO: shouldn't we add the function prototype to the header file? // Sets the rows of constraints matrix G for the energy-conserving EQP rule. void ComputeRowsOfG(const IntegrationRule *ir, hydrodynamics::QuadratureData const& quad_data, - Vector const& w_e, Vector const& v_e, + Vector const& v_j_e, Vector const& w_j_e, Vector const& v_i_e, FiniteElement const& test_fe, FiniteElement const& trial_fe, const int zone_id, @@ -357,26 +357,28 @@ void ComputeRowsOfG(const IntegrationRule *ir, if (equationV) { MFEM_VERIFY(rv.Size() == nqp, ""); + MFEM_VERIFY(v_j_e.Size() == h1dofs_cnt*dim, ""); } if (equationE) { MFEM_VERIFY(re.Size() == nqp, ""); - MFEM_VERIFY(v_e.Size() == h1dofs_cnt*dim, ""); - MFEM_VERIFY(w_e.Size() == l2dofs_cnt, ""); + MFEM_VERIFY(v_i_e.Size() == h1dofs_cnt*dim, ""); + MFEM_VERIFY(w_j_e.Size() == l2dofs_cnt, ""); } - DenseMatrix grad_vshape(h1dofs_cnt, dim); // grad of velocity basis vector + DenseMatrix grad_vshape(h1dofs_cnt, dim); // grad of velocity shape function DenseMatrix loc_force(h1dofs_cnt, dim); Vector Vloc_force(loc_force.Data(), h1dofs_cnt * dim); + Vector eshape(l2dofs_cnt), unitE(l2dofs_cnt); + unitE = 1.0; + for (int q = 0; q < nqp; q++) { const IntegrationPoint &ip = ir->IntPoint(q); // Form stress:grad_vshape at the current point. - // TODO: how do we set which velocity basis vector is used for the - // gradient calculation? test_fe.CalcDShape(ip, grad_vshape); for (int i = 0; i < h1dofs_cnt; i++) { @@ -394,25 +396,18 @@ void ComputeRowsOfG(const IntegrationRule *ir, // NOTE: UpdateQuadratureData includes ip.weight as a factor in quad_data.stressJinvT, // set by LagrangianHydroOperator::UpdateQuadratureData. loc_force *= 1.0 / ip.weight; // Divide by exact quadrature weight + + trial_fe.CalcShape(ip, eshape); // Energy shape function if (equationV) { - // set rv = sum(Vloc_force) - rv[q] = 0.0; - for (int i = 0; i < h1dofs_cnt * dim; i++) - { - rv[q] += Vloc_force(i); - } + // Inner product, on this element, with the jth V basis vector. + rv[q] = (v_j_e * Vloc_force) * (eshape * unitE); } if (equationE) { - // set re = F^T * v; F = integral of (stress : grad w_i)phi_j - re[q] = 0.0; - for (int i = 0; i < l2dofs_cnt; i++) - { - re[q] += w_e(i); - } - re[q] *= v_e * Vloc_force; + // Inner product, on this element, with the jth E basis vector. + re[q] = (v_i_e * Vloc_force) * (eshape * w_j_e); } } // q -- quadrature point } @@ -747,8 +742,7 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, CAROM::Matrix Wv(H1size, NBv, true); CAROM::Matrix We(L2size, NBe, true); - // jth Wv columnn holds the jth velocity basis vector - // TODO: do we need Wv anymore? it's not used anywhere + // jth Wv columnn holds the jth velocity ROM basis vector for (int j=0; j const& w_el = ir0->GetWeights(); MFEM_VERIFY(w_el.Size() == nqe, ""); + ParGridFunction gf2H1(gfH1); + for (int i=0; iGetWeights(). - // In the above, w_j is the jth velocity basis vector. + // In the above, e_j is the jth velocity basis vector, 1E is the + // identity vector in the energy space. // Energy equation. // For 0 <= j < NBe, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, // entry G(estart + j + (i*NBe), (e*nqe) + m) is the coefficient of // - // (stress(v_i,e_i,x_i) : grad v_i) w_j + // e_j^T * F(v_i,e_i,x_i)^T * v_i // // at point m of element e with respect to the integration rule weight at // that point, where the "exact" quadrature solution is ir0->GetWeights(). - // In the above, w_j is the jth energy basis vector. + // In the above, e_j is the jth energy basis vector, v_i is the ith + // velocity snapshot. // Set the contraints for velocity and energy at the same time, then // add the rest for velocity or energy, depending on which variable has @@ -839,7 +837,7 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, for (int j=0; jnumRows()); for (int k = 0; k < basisE->numRows(); ++k) @@ -847,25 +845,25 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, Ediff[k] = We(k, j); } gfL2.SetFromTrueDofs(Ediff); // jth energy basis vector - gfH1 = S_v; // ith velocity snapshot + gf2H1 = S_v; // ith velocity snapshot for (int e=0; eGetQuadData(), - w_j_e, v_i_e, + v_j_e, w_j_e, v_i_e, *input.H1FESpace->GetFE(e), *input.L2FESpace->GetFE(e), e, true, true, rv, re); @@ -882,18 +880,18 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, for (int j=NBmin; jGetQuadData(), - w_j_e, v_i_e, + v_j_e, w_j_e, v_i_e, *input.H1FESpace->GetFE(e), *input.L2FESpace->GetFE(e), e, true, false, rv, re); @@ -911,22 +909,23 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, Vector Ediff(basisE->numRows()); for (int k = 0; k < basisE->numRows(); ++k) { - Ediff[k] = We(k, j); // jth energy basis vector + Ediff[k] = We(k, j); } - gfL2.SetFromTrueDofs(Ediff); - gfH1 = S_v; // ith velocity snapshot + gfL2.SetFromTrueDofs(Ediff); // jth energy basis vector + gf2H1 = S_v; // ith velocity snapshot for (int e=0; eGetQuadData(), - w_j_e, v_i_e, + v_j_e, w_j_e, v_i_e, *input.H1FESpace->GetFE(e), *input.L2FESpace->GetFE(e), e, false, true, rv, re); @@ -2780,6 +2779,7 @@ ROM_Operator::ROM_Operator(ROM_Options const& input, ROM_Basis *b, fy.SetSize(input.FOMoper->Height()); } + // TODO: remove this and the useReducedM flag? if (useReducedM) { ComputeReducedMv(); @@ -2788,15 +2788,12 @@ ROM_Operator::ROM_Operator(ROM_Options const& input, ROM_Basis *b, if (hyperreduce && hyperreductionSamplingType == eqp) { - // basic EQP - // TODO: are reduced mass matrices needed for EQP, or just W matrices? - - // TODO: why compute the reduced mass matrices, when for the basic EQP - // case we've already multiplied my M^{-1} and basis^T when deriving - // the reduced quadrature rules? + // basic EQP - // ComputeReducedMv(); - // ComputeReducedMe(); + // Computes the product Phi_v^T * Mv^{-1} (similarly for energy) + // used in forming the RHS force vectors. + ComputeReducedMv(); + ComputeReducedMe(); ReadSolutionNNLS(input, "run/nnlsV", eqpI, eqpW); ReadSolutionNNLS(input, "run/nnlsE", eqpI_E, eqpW_E); @@ -2804,8 +2801,11 @@ ROM_Operator::ROM_Operator(ROM_Options const& input, ROM_Basis *b, else if (hyperreduce && hyperreductionSamplingType == eqp_energy) { // energy-conserving EQP - // NOTE: will not compute reduced M matrices; instead, will normalize - // the basis matrices so that the reduced M matrices are identity. + + // TODO: use them to compute Phi_v^T * Mv * Phi_v (similarly for energy) + // needed to form the RHS vectors to time march. + ComputeReducedMv(); + ComputeReducedMe(); ReadSolutionNNLS(input, "run/nnlsEC", eqpI, eqpW); } @@ -3177,6 +3177,10 @@ void ROM_Operator::ComputeReducedMv() (*Wmat)(i,j) = Mvj[i]; } } + else if (hyperreduce && hyperreductionSamplingType == eqp_energy) + { + // TODO: form the inverse of the reduced mass matrix + } else if (!hyperreduce) { MFEM_ABORT("TODO"); @@ -3224,6 +3228,10 @@ void ROM_Operator::ComputeReducedMe() (*Wmat_E)(i,j) = Mej[i]; } } + else if (hyperreduce && hyperreductionSamplingType == eqp_energy) + { + // TODO: form the inverse of the reduced mass matrix + } else if (!hyperreduce) { MFEM_ABORT("TODO"); @@ -3467,7 +3475,6 @@ void ROM_Operator::UndoInducedGramSchmidt(const int var, Vector &S, bool keep_da void ROM_Operator::ApplyHyperreduction(Vector &S) { - // TODO: should set GramSchmidt = true when using energy-conserving EQP. if (useGramSchmidt && !sns1) { InducedGramSchmidt(1, S); // velocity @@ -4463,7 +4470,7 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res, const int nqe = ir->GetWeights().Size(); - DenseMatrix vshape, loc_force; + DenseMatrix grad_vshape, loc_force; Array vdofs; res = 0.0; @@ -4558,11 +4565,11 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res, const int l2dofs_cnt = trial_fe->GetDof(); - vshape.SetSize(h1dofs_cnt, dim); + grad_vshape.SetSize(h1dofs_cnt, dim); loc_force.SetSize(h1dofs_cnt, dim); - // Form stress:grad_shape at the current point. - test_fe->CalcDShape(ip, vshape); + // Form stress:grad_vshape at the current point. + test_fe->CalcDShape(ip, grad_vshape); for (int k = 0; k < h1dofs_cnt; k++) { for (int vd = 0; vd < dim; vd++) // Velocity components. @@ -4571,7 +4578,7 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res, for (int gd = 0; gd < dim; gd++) // Gradient components. { loc_force(k, vd) += - quad_data.stressJinvT(vd)(e*nqe + qpi, gd) * vshape(k,gd); + quad_data.stressJinvT(vd)(e*nqe + qpi, gd) * grad_vshape(k,gd); } } } @@ -4580,6 +4587,10 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res, Vector Vloc_force(loc_force.Data(), loc_force.NumRows() * loc_force.NumCols()); Vector v_e(h1dofs_cnt * dim); + Vector eshape(l2dofs_cnt), unitE(l2dofs_cnt); + trial_fe->CalcShape(ip, eshape); + unitE = 1.0; + if (energy_conserve) { // energy-conserving EQP @@ -4587,30 +4598,24 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res, { // TODO: v_e should be the jth velocity basis vector DOFs that // correspond to the current element. - // Recall that the basis vectors have been renormalized by the - // induced Gram-Schmidt process, so that reduced Mv = I. - // I think the below is not right. basis->GetBasisVectorV(false, j, v_e); - res[j] += v_e * Vloc_force; + // Inner product, on this element, with the j-th V basis vector. + res[j] += (v_e * Vloc_force) * (eshape * unitE); } } else { // basic EQP - Vector shape(l2dofs_cnt), unitE(l2dofs_cnt); - trial_fe->CalcShape(ip, shape); - unitE = 1.0; - const int eos = elemIndex * nvdof; for (int j = 0; j < rdim; ++j) { for (int k = 0; k < nvdof; ++k) v_e[k] = W_elems(eos + k, j); - // Compute the inner product, on this element, with the j-th W vector. - // TODO: why compute shape * unitE just to get 1.0 back? - res[j] += (v_e * Vloc_force) * (shape * unitE); + // Inner product, on this element, with the j-th W vector. + // W is the product of Phi_v^T * Mv^{-1} + res[j] += (v_e * Vloc_force) * (eshape * unitE); } } } // Loop (i) over EQP points @@ -4629,7 +4634,7 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res, const int nqe = ir->GetWeights().Size(); - DenseMatrix vshape, loc_force; + DenseMatrix grad_vshape, loc_force; Array vdofs, edofs; Vector v_e; @@ -4741,11 +4746,11 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res, const int l2dofs_cnt = trial_fe->GetDof(); - vshape.SetSize(h1dofs_cnt, dim); + grad_vshape.SetSize(h1dofs_cnt, dim); loc_force.SetSize(h1dofs_cnt, dim); - // Form stress:grad_shape at the current point. - test_fe->CalcDShape(ip, vshape); + // Form stress:grad_vshape at the current point. + test_fe->CalcDShape(ip, grad_vshape); for (int k = 0; k < h1dofs_cnt; k++) { for (int vd = 0; vd < dim; vd++) // Velocity components. @@ -4754,7 +4759,7 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res, for (int gd = 0; gd < dim; gd++) // Gradient components. { loc_force(k, vd) += - quad_data.stressJinvT(vd)(e*nqe + qpi, gd) * vshape(k,gd); + quad_data.stressJinvT(vd)(e*nqe + qpi, gd) * grad_vshape(k,gd); } } } @@ -4763,31 +4768,33 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res, Vector Vloc_force(loc_force.Data(), loc_force.NumRows() * loc_force.NumCols()); - Vector shape(l2dofs_cnt); - trial_fe->CalcShape(ip, shape); + Vector eshape(l2dofs_cnt), w_e(nedof); + trial_fe->CalcShape(ip, eshape); if (energy_conserve) { // energy-conserving EQP for (int j = 0; j < rdim; j++) { - for (int k = 0; k < l2dofs_cnt; k++) res[j] += shape(k); + // TODO: w_e should be the jth energy basis vector DOFs that + // correspond to the current element. - res[j] *= v_e * Vloc_force; + // Inner product, on this element, with the jth E basis vector. + res[j] += (v_e * Vloc_force) * (eshape * w_e); } } else { // basic EQP - Vector w_e(nedof); const int eos = elemIndex * nedof; for (int j=0; jForceIntegratorEQP_FOM(rhs, rom_op->getSamplingType() == eqp_energy); } else @@ -410,8 +408,6 @@ void LagrangianHydroOperator::SolveEnergy(const Vector &S, const Vector &v, timer.sw_force.Start(); if (use_eqp) { - // TODO: this function doesn't have access to the hyperreduction - // sampling type; how should I resolve this? rom_op->ForceIntegratorEQP_E_FOM(v, e_rhs, rom_op->getSamplingType() == eqp_energy); } else From 36b2232e21130da695ec0e46ade71903c9577e70 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Tue, 11 Jul 2023 17:11:06 -0700 Subject: [PATCH 09/39] form inverse reduced mass matrices --- rom/laghos_rom.cpp | 50 ++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 48 insertions(+), 2 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 6453a56e..0c830062 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -3177,9 +3177,34 @@ void ROM_Operator::ComputeReducedMv() (*Wmat)(i,j) = Mvj[i]; } } + // TODO: do I need to enforce MPI rank == 0? else if (hyperreduce && hyperreductionSamplingType == eqp_energy) { - // TODO: form the inverse of the reduced mass matrix + // Form inverse of reduced Mv + invMvROM.SetSize(nv); + + const int size_H1 = basis->SolutionSizeH1FOM(); + const int tsize_H1 = H1spaceFOM->GetTrueVSize(); + + ParGridFunction gf(H1spaceFOM), gf2(H1spaceFOM); + Vector vj(tsize_H1), vi(tsize_H1); + Vector Mvj(size_H1); + + for (int j = 0; j < nv; ++j) + { + basis->GetBasisVectorV(false, j, vj); + gf.SetFromTrueDofs(vj); + operFOM->MultMv(gf, Mvj); + + for (int i = 0; i < nv; ++i) + { + basis->GetBasisVectorV(false, i, vi); + gf2.SetFromTrueDofs(vi); + + invMvROM(i,j) = gf2 * Mvj; + } + } + invMvROM.Invert(); } else if (!hyperreduce) { @@ -3228,9 +3253,30 @@ void ROM_Operator::ComputeReducedMe() (*Wmat_E)(i,j) = Mej[i]; } } + // TODO: do I need to enforce MPI rank == 0? else if (hyperreduce && hyperreductionSamplingType == eqp_energy) { - // TODO: form the inverse of the reduced mass matrix + // Form inverse of reduced Me + invMeROM.SetSize(ne); + + const int size_L2 = basis->SolutionSizeL2FOM(); + + Vector ej(size_L2), ei(size_L2); + Vector Mej(size_L2); + + for (int j = 0; j < ne; ++j) + { + basis->GetBasisVectorE(false, j, ej); + operFOM->MultMe(ej, Mej); + + for (int i = 0; i < ne; ++i) + { + basis->GetBasisVectorE(false, i, ei); + + invMvROM(i,j) = ei * Mej; + } + } + invMvROM.Invert(); } else if (!hyperreduce) { From 684a03985962323d534466feb53bb1e43319b097 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Wed, 12 Jul 2023 13:37:16 -0700 Subject: [PATCH 10/39] form RHSs for timemarching --- rom/laghos_rom.cpp | 64 +++++++++++++++++++++++++++++++++++----------- 1 file changed, 49 insertions(+), 15 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 0d1a6ebe..0ff5579b 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -339,6 +339,7 @@ void ComputeElementRowOfG_E(const IntegrationRule *ir, } // TODO: shouldn't we add the function prototype to the header file? + // Sets the rows of constraints matrix G for the energy-conserving EQP rule. void ComputeRowsOfG(const IntegrationRule *ir, hydrodynamics::QuadratureData const& quad_data, @@ -2802,8 +2803,8 @@ ROM_Operator::ROM_Operator(ROM_Options const& input, ROM_Basis *b, { // energy-conserving EQP - // TODO: use them to compute Phi_v^T * Mv * Phi_v (similarly for energy) - // needed to form the RHS vectors to time march. + // Computes Phi_v^T * Mv * Phi_v (similarly for energy) needed to + // form the RHS vectors to time march. ComputeReducedMv(); ComputeReducedMe(); @@ -3185,6 +3186,8 @@ void ROM_Operator::ComputeReducedMv() const int size_H1 = basis->SolutionSizeH1FOM(); const int tsize_H1 = H1spaceFOM->GetTrueVSize(); + + Wmat = new CAROM::Matrix(size_H1, nv, true); ParGridFunction gf(H1spaceFOM), gf2(H1spaceFOM); Vector vj(tsize_H1), vi(tsize_H1); @@ -3194,6 +3197,10 @@ void ROM_Operator::ComputeReducedMv() { basis->GetBasisVectorV(false, j, vj); gf.SetFromTrueDofs(vj); + + // store jth V basis vector in Wmat(:,j) + for (int i=0; iMultMv(gf, Mvj); for (int i = 0; i < nv; ++i) @@ -3261,12 +3268,18 @@ void ROM_Operator::ComputeReducedMe() const int size_L2 = basis->SolutionSizeL2FOM(); + Wmat_E = new CAROM::Matrix(size_L2, ne, true); + Vector ej(size_L2), ei(size_L2); Vector Mej(size_L2); for (int j = 0; j < ne; ++j) { basis->GetBasisVectorE(false, j, ej); + + // store jth E basis vector in Wmat_E(:,j) + for (int i=0; iMultMe(ej, Mej); for (int i = 0; i < ne; ++i) @@ -4647,35 +4660,45 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res, trial_fe->CalcShape(ip, eshape); unitE = 1.0; + const int eos = elemIndex * nvdof; + if (energy_conserve) { // energy-conserving EQP for (int j = 0; j < rdim; ++j) { - // TODO: v_e should be the jth velocity basis vector DOFs that - // correspond to the current element. - // I think the below is not right. - basis->GetBasisVectorV(false, j, v_e); + // v_e: jth V basis vector's DOFs on this element + for (int k = 0; k < nvdof; ++k) v_e[k] = W_elems(eos + k, j); - // Inner product, on this element, with the j-th V basis vector. + // Inner product, on this element, with the jth V basis vector. res[j] += (v_e * Vloc_force) * (eshape * unitE); } } else { // basic EQP - const int eos = elemIndex * nvdof; for (int j = 0; j < rdim; ++j) { + // v_e is the product Phi_v^T * Mv^{-1} for (int k = 0; k < nvdof; ++k) v_e[k] = W_elems(eos + k, j); - // Inner product, on this element, with the j-th W vector. - // W is the product of Phi_v^T * Mv^{-1} + // Inner product, on this element, with the jth W vector. res[j] += (v_e * Vloc_force) * (eshape * unitE); } } } // Loop (i) over EQP points + if (energy_conserve) + { + // Multiply by the reduced Mv inverse + + // TODO: does that copy the data of res to res_tmp, or is it just + // a reference? + Vector res_tmp = res; + + invMvROM.Mult(res_tmp, res); + } + MPI_Allreduce(MPI_IN_PLACE, res.GetData(), res.Size(), MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); } @@ -4827,13 +4850,15 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res, Vector eshape(l2dofs_cnt), w_e(nedof); trial_fe->CalcShape(ip, eshape); + const int eos = elemIndex * nedof; + if (energy_conserve) { // energy-conserving EQP for (int j = 0; j < rdim; j++) { - // TODO: w_e should be the jth energy basis vector DOFs that - // correspond to the current element. + // w_e: jth E basis vector's DOFs on this element + for (int k=0; k Date: Fri, 14 Jul 2023 15:22:12 -0700 Subject: [PATCH 11/39] bug fixes & minor changes --- rom/laghos.cpp | 28 ++++++++++++++-------------- rom/laghos_rom.cpp | 20 ++++++++------------ 2 files changed, 22 insertions(+), 26 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 2a100947..14a17ec5 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -395,13 +395,13 @@ int main(int argc, char *argv[]) romOptions.hyperreductionSamplingType = getHyperreductionSamplingType(hyperreductionSamplingType); romOptions.use_sample_mesh = romOptions.hyperreduce && (romOptions.hyperreductionSamplingType != eqp - || romOptions.hyperreductionSamplingType != eqp_energy); + && romOptions.hyperreductionSamplingType != eqp_energy); MFEM_VERIFY(!(romOptions.SNS) || (romOptions.hyperreductionSamplingType != eqp && romOptions.hyperreductionSamplingType != eqp_energy), "Using SNS with EQP is prohibited"); - if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) + if (romOptions.hyperreductionSamplingType == eqp_energy) romOptions.GramSchmidt = false; romOptions.spaceTimeMethod = getSpaceTimeMethod(spaceTimeMethod); @@ -1004,14 +1004,14 @@ int main(int argc, char *argv[]) LagrangianHydroOperator* oper = NULL; if (fom_data) { - const bool noMassSolve = rom_online && (romOptions.hyperreductionSamplingType == eqp); + const bool noMassSolve = rom_online && (romOptions.hyperreductionSamplingType == eqp || romOptions.hyperreductionSamplingType == eqp_energy); oper = new LagrangianHydroOperator(S->Size(), *H1FESpace, *L2FESpace, ess_tdofs, *rho, source, cfl, mat_gf_coeff, visc, vort, p_assembly, cg_tol, cg_max_iter, ftz_tol, H1FEC.GetBasisType(), noMassSolve, noMassSolve, - rom_online && (romOptions.hyperreductionSamplingType == eqp)); + rom_online && (romOptions.hyperreductionSamplingType == eqp || romOptions.hyperreductionSamplingType == eqp_energy)); } socketstream* vis_rho = NULL; @@ -1581,11 +1581,11 @@ int main(int argc, char *argv[]) romOper[0]->ApplyHyperreduction(romS); } - // TODO: do we want that for the energy-conserving EQP? - if (rom_online && romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) - { - romOper[0]->ApplyHyperreduction(romS); - } + // // TODO: do we want that for the energy-conserving EQP? + // if (rom_online && romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) + // { + // romOper[0]->ApplyHyperreduction(romS); + // } double windowEndpoint = 0.0; double windowOverlapMidpoint = 0.0; @@ -2058,10 +2058,10 @@ int main(int argc, char *argv[]) romOper[romOptions.window]->ApplyHyperreduction(romS); } - if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) - { - romOper[romOptions.window]->ApplyHyperreduction(romS); - } + //if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) + //{ + // romOper[romOptions.window]->ApplyHyperreduction(romS); + //} if (problem == 7 && romOptions.indicatorType == penetrationDistance) { @@ -2109,7 +2109,7 @@ int main(int argc, char *argv[]) << sqrt(tot_norm) << endl; } - if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) + if (romOptions.hyperreductionSamplingType == eqp_energy) { double energy_total, energy_diff; energy_total = oper->InternalEnergy(*e_gf) + diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 0ff5579b..b385c035 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -1119,16 +1119,10 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) // Form the energy identity and replace the last basis vector by it. Vector unitE(tL2size); - // TODO: my understanding is that the "=" opetaror is overloaded, so - // that, by setting unitE = 1.0, the vector unitE contains the coeffs - // that result in the summed function being 1.0 on the grid points. - // Is that correct? unitE = 1.0; for (int i = 0; i < tL2size; i++) { - // TODO: I believe I can access the last column in this way, - // without having to use pointer arithmetic, right? (*tBasisE)(i, cutoff[2]-1) = unitE[i]; } @@ -2808,7 +2802,12 @@ ROM_Operator::ROM_Operator(ROM_Options const& input, ROM_Basis *b, ComputeReducedMv(); ComputeReducedMe(); + // Read the same data twice, because different variables are used + // when solving the velocity and energy problems (due to the way this + // is done for basic EQP). In this way, minimal changes to the code + // are needed. ReadSolutionNNLS(input, "run/nnlsEC", eqpI, eqpW); + ReadSolutionNNLS(input, "run/nnlsEC", eqpI_E, eqpW_E); } } @@ -3286,10 +3285,10 @@ void ROM_Operator::ComputeReducedMe() { basis->GetBasisVectorE(false, i, ei); - invMvROM(i,j) = ei * Mej; + invMeROM(i,j) = ei * Mej; } } - invMvROM.Invert(); + invMeROM.Invert(); } else if (!hyperreduce) { @@ -4691,9 +4690,6 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res, if (energy_conserve) { // Multiply by the reduced Mv inverse - - // TODO: does that copy the data of res to res_tmp, or is it just - // a reference? Vector res_tmp = res; invMvROM.Mult(res_tmp, res); @@ -4907,7 +4903,7 @@ void ROM_Operator::ForceIntegratorEQP_E_FOM(Vector const& v, Vector & rhs, { Vector res(basis->GetDimE()); - ForceIntegratorEQP_E(v, res); + ForceIntegratorEQP_E(v, res, energy_conserve); basis->LiftROMtoFOM_dEdt(res, rhs); } From 679c952196d20570a05d6efce50b751408ada642 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Tue, 25 Jul 2023 09:20:24 -0700 Subject: [PATCH 12/39] fix energy identity mistake --- rom/laghos.cpp | 12 ++++++++---- rom/laghos_rom.cpp | 42 +++++++++++++++++++++++++++++++----------- 2 files changed, 39 insertions(+), 15 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 14a17ec5..af66ef6e 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -2118,9 +2118,9 @@ int main(int argc, char *argv[]) if (mpi.Root()) { - cout << "\tE_tot = " << setprecision(14) + cout << "\tE_tot = " << scientific << setprecision(5) << energy_total - << ",\tE_diff = " << setprecision(14) + << ",\tE_diff = " << scientific << setprecision(5) << energy_diff << endl; } } @@ -2502,8 +2502,12 @@ int main(int argc, char *argv[]) if (mpi.Root()) { cout << endl; - cout << "Energy diff: " << scientific << setprecision(2) - << fabs(energy_init - energy_final) << endl; + cout << "Initial energy: " << scientific << setprecision(5) + << energy_init << endl; + cout << "Energy diff: " << scientific << setprecision(5) + << energy_final - energy_init << endl; + cout << "Rel. energy diff: " << scientific << setprecision(5) + << (energy_final - energy_init) / energy_init << endl; } PrintParGridFunction(myid, testing_parameter_outputPath + "/x_gf" + romOptions.basisIdentifier, x_gf); diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index b385c035..1b09183a 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -988,16 +988,6 @@ void ROM_Sampler::SetupEQP_Force(const CAROM::Matrix* snapX, void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) { - if (input.hyperreductionSamplingType == eqp_energy) - { - // For the energy-conserving EQP case, we increase the energy basis - // dimension by 1 to accomodate for the addition of the energy - // identity. - // The change is done here, before the window basis parameters are - // written to their files, and is also carried to the caller's scope. - cutoff[2] += 1; - } - if (writeSnapshots) { if (!useXV) generator_X->writeSnapshot(); @@ -1105,6 +1095,22 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) if (input.hyperreductionSamplingType == eqp_energy) { + if (rank == 0) + { + // For the energy-conserving EQP case, we increase the energy basis + // dimension by 1 to accomodate for the addition of the energy + // identity. + // The change is done here, before the window basis parameters are + // written to their files, and is also carried to the caller's scope. + + // TODO: if cutoff[2] == windowNumSamples then increasing this by + // 1 will lead to a function call trying to get more columns than + // there are available in basisE. In this case, we'd need to take + // all columns of basisE and then append one. + // And then do the same during the online stage. + cutoff[2] += 1; + } + const CAROM::Matrix *basisV = generator_V->getSpatialBasis(); const CAROM::Matrix *basisE = generator_E->getSpatialBasis(); @@ -1118,7 +1124,6 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) // Form the energy identity and replace the last basis vector by it. Vector unitE(tL2size); - unitE = 1.0; for (int i = 0; i < tL2size; i++) @@ -1126,6 +1131,8 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) (*tBasisE)(i, cutoff[2]-1) = unitE[i]; } + tBasisE->orthogonalize(); + SetupEQP_Force(generator_X->getSnapshotMatrix(), generator_V->getSnapshotMatrix(), generator_E->getSnapshotMatrix(), @@ -1282,6 +1289,19 @@ ROM_Basis::ROM_Basis(ROM_Options const& input, MPI_Comm comm_, const double sFac mfL2.SetSize(tL2size); ReadSolutionBases(input.window); + + if (hyperreduce && hyperreductionSamplingType == eqp_energy) + { + // Form the energy identity and replace the last E-basis vector by it. + Vector unitE(tL2size); + unitE = 1.0; + for (int i = 0; i < tL2size; i++) + { + (*basisE)(i, rdime-1) = unitE[i]; + } + basisE->orthogonalize(); + } + if (spaceTime) { ReadTemporalBases(input.window); From c2c7fd5ed7b6079f956fc0960917a8bdccf4f8f9 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Tue, 25 Jul 2023 18:12:03 -0700 Subject: [PATCH 13/39] fix variables offset --- rom/laghos_rom.cpp | 131 +++++++++++++++++++++++++-------------------- rom/laghos_rom.hpp | 5 +- 2 files changed, 77 insertions(+), 59 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 1b09183a..0eeab94b 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -57,9 +57,10 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p SetStateVariables(S); SetStateVariableRates(dt); - const bool sampleX = generator_X->isNextSample(t); - - Vector dSdt; + const bool sampleX = (hyperreductionSamplingType == eqp_energy) ? true : + generator_X->isNextSample(t); + + Vector dSdt; if (!sns && rhsBasis) { dSdt.SetSize(S.Size()); @@ -98,7 +99,8 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p } } - const bool sampleV = generator_V->isNextSample(t); + const bool sampleV = (hyperreductionSamplingType == eqp_energy) ? true : + generator_V->isNextSample(t); //TODO: use this, plus generator_Fv->computeNextSampleTime? So far, it seems sampleV == true on every step. //const bool sampleFv = generator_Fv->isNextSample(t); @@ -152,7 +154,8 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p } } - const bool sampleE = generator_E->isNextSample(t); + const bool sampleE = (hyperreductionSamplingType == eqp_energy) ? true : + generator_E->isNextSample(t); if (sampleE) { @@ -672,7 +675,7 @@ void ROM_Sampler::SetupEQP_Force_Eq(const CAROM::Matrix* snapX, } // j } // i - CAROM::Vector w(ne * nqe, true); + CAROM::Vector w(ne * nqe, true); for (int i=0; iResetQuadratureData(); // Velocity equation. - // For 0 <= j < NBv, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, + // For 0 <= j < NBv, 0 <= i < nsnap, 0 <= e < ne, 0 <= m < nqe, // entry G(j + (i*NBv), (e*nqe) + m) is the coefficient of // // e_j^T * F(v_i,e_i,x_i) * 1E @@ -821,7 +816,7 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, // identity vector in the energy space. // Energy equation. - // For 0 <= j < NBe, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe, + // For 0 <= j < NBe, 0 <= i < nsnap, 0 <= e < ne, 0 <= m < nqe, // entry G(estart + j + (i*NBe), (e*nqe) + m) is the coefficient of // // e_j^T * F(v_i,e_i,x_i)^T * v_i @@ -1097,17 +1092,10 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) { if (rank == 0) { - // For the energy-conserving EQP case, we increase the energy basis + // For the energy-conserving EQP case, increase the energy basis // dimension by 1 to accomodate for the addition of the energy // identity. - // The change is done here, before the window basis parameters are - // written to their files, and is also carried to the caller's scope. - // TODO: if cutoff[2] == windowNumSamples then increasing this by - // 1 will lead to a function call trying to get more columns than - // there are available in basisE. In this case, we'd need to take - // all columns of basisE and then append one. - // And then do the same during the online stage. cutoff[2] += 1; } @@ -1120,17 +1108,21 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) // For the energy basis, the cutoff parameter has already been // increased by 1 to accomodate the energy identity. CAROM::Matrix *tBasisV = basisV->getFirstNColumns(cutoff[1]); - CAROM::Matrix *tBasisE = basisE->getFirstNColumns(cutoff[2]); - - // Form the energy identity and replace the last basis vector by it. - Vector unitE(tL2size); - unitE = 1.0; + + CAROM::Matrix *tBasisE = new CAROM::Matrix(tL2size, cutoff[2], false); + // Get the first cutoff[2]-1 columns of basisE for (int i = 0; i < tL2size; i++) { - (*tBasisE)(i, cutoff[2]-1) = unitE[i]; + for (int j = 0; j < cutoff[2] - 1; j++) + (*tBasisE)(i, j) = (*basisE)(i, j); } + // Form the energy identity and include it as the last basis vector. + Vector unitE(tL2size); + unitE = 1.0; + for (int i = 0; i < tL2size; i++) + (*tBasisE)(i, cutoff[2]-1) = unitE[i]; tBasisE->orthogonalize(); SetupEQP_Force(generator_X->getSnapshotMatrix(), @@ -1288,17 +1280,18 @@ ROM_Basis::ROM_Basis(ROM_Options const& input, MPI_Comm comm_, const double sFac mfH1.SetSize(tH1size); mfL2.SetSize(tL2size); + if (hyperreduce && hyperreductionSamplingType == eqp_energy) + basisE = new CAROM::Matrix(tL2size, rdime, true); + ReadSolutionBases(input.window); if (hyperreduce && hyperreductionSamplingType == eqp_energy) { - // Form the energy identity and replace the last E-basis vector by it. + // Form the energy identity and include it as the last basis vector. Vector unitE(tL2size); unitE = 1.0; for (int i = 0; i < tL2size; i++) - { (*basisE)(i, rdime-1) = unitE[i]; - } basisE->orthogonalize(); } @@ -2222,7 +2215,29 @@ void ROM_Basis::ReadSolutionBases(const int window) if (!useVX) basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + std::to_string(window) + basisIdentifier, tH1size, rdimv); - basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + std::to_string(window) + basisIdentifier, tL2size, rdime); + // In the energy-conserving EQP case we read the first rdime-1 basis + // vectors, since the rdime parameter has been increased by 1 to + // accommodate the addition of the energy identity. + if (hyperreduce && hyperreductionSamplingType == eqp_energy) + { + int tmp_rdime = rdime - 1; + CAROM::Matrix *tmp_basisE = 0; + + tmp_basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + + std::to_string(window) + basisIdentifier, tL2size, tmp_rdime); + + for (int i = 0; i < tL2size; i++) + { + for (int j = 0; j < tmp_rdime; j++) + (*basisE)(i, j) = (*tmp_basisE)(i, j); + } + } + else + { + basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + + std::to_string(window) + basisIdentifier, tL2size, rdime); + } + if (useXV) basisX = basisV; diff --git a/rom/laghos_rom.hpp b/rom/laghos_rom.hpp index f9b6fc4e..87fe723a 100644 --- a/rom/laghos_rom.hpp +++ b/rom/laghos_rom.hpp @@ -481,7 +481,8 @@ class ROM_Sampler energyFraction_X(input.energyFraction_X), sns(input.SNS), lhoper(input.FOMoper), writeSnapshots(input.parameterID >= 0), parameterID(input.parameterID), basename(*input.basename), Voffset(!input.useXV && !input.useVX && !input.mergeXV), useXV(input.useXV), useVX(input.useVX), VTos(input.VTos), spaceTime(input.spaceTimeMethod != no_space_time), - rhsBasis(input.hyperreductionSamplingType != eqp && input.hyperreductionSamplingType != eqp_energy) + rhsBasis(input.hyperreductionSamplingType != eqp && input.hyperreductionSamplingType != eqp_energy), + hyperreductionSamplingType(input.hyperreductionSamplingType) { const int window = input.window; @@ -685,6 +686,8 @@ class ROM_Sampler const bool spaceTime; + const bool hyperreductionSamplingType; + hydrodynamics::LagrangianHydroOperator *lhoper; int VTos = 0; // Velocity temporal index offset, used for V and Fe. This fixes the issue that V and Fe are not sampled at t=0, since they are initially zero. This is valid for the Sedov test but not in general when the initial velocity is nonzero. From 4ade8619bb1867f3885cca76c9a3fadcd1f047fc Mon Sep 17 00:00:00 2001 From: Chris Val Date: Mon, 31 Jul 2023 15:36:31 -0700 Subject: [PATCH 14/39] fix restore issue --- rom/laghos_rom.cpp | 36 ++++++++++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 0eeab94b..ed1438f8 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -934,6 +934,9 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, } // j -- basis vector } // i -- snapshot + // Rescale every column (NNLS equation) by its max absolute value + //Gt.rescale_cols_max(); + CAROM::Vector w(ne * nqe, true); for (int i=0; i tmpmax) + { + tmpmax = fabs(Gt(ii, jj)); + } + } + + if (tmpmax > 1.0e-14) + { + for (int ii = 0; ii < nGtrows; ii++) + { + Gt(ii, jj) /= tmpmax; + } + } + } + } + // TODO: input these NNLS parameters? double tolNNLS = 1.0e-14; @@ -1093,7 +1121,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) if (rank == 0) { // For the energy-conserving EQP case, increase the energy basis - // dimension by 1 to accomodate for the addition of the energy + // dimension by 1 to accomodate the addition of the energy // identity. cutoff[2] += 1; @@ -1280,12 +1308,12 @@ ROM_Basis::ROM_Basis(ROM_Options const& input, MPI_Comm comm_, const double sFac mfH1.SetSize(tH1size); mfL2.SetSize(tL2size); - if (hyperreduce && hyperreductionSamplingType == eqp_energy) + if (hyperreductionSamplingType == eqp_energy) basisE = new CAROM::Matrix(tL2size, rdime, true); ReadSolutionBases(input.window); - if (hyperreduce && hyperreductionSamplingType == eqp_energy) + if (hyperreductionSamplingType == eqp_energy) { // Form the energy identity and include it as the last basis vector. Vector unitE(tL2size); @@ -2218,7 +2246,7 @@ void ROM_Basis::ReadSolutionBases(const int window) // In the energy-conserving EQP case we read the first rdime-1 basis // vectors, since the rdime parameter has been increased by 1 to // accommodate the addition of the energy identity. - if (hyperreduce && hyperreductionSamplingType == eqp_energy) + if (hyperreductionSamplingType == eqp_energy) { int tmp_rdime = rdime - 1; CAROM::Matrix *tmp_basisE = 0; From 834de51c4ca0d09a5e9e9654f8c58f08055a3caf Mon Sep 17 00:00:00 2001 From: Chris Val Date: Fri, 28 Jul 2023 16:03:20 -0700 Subject: [PATCH 15/39] add NNLS tolerance input option --- rom/laghos.cpp | 7 +++++-- rom/laghos_rom.cpp | 12 +++--------- rom/laghos_rom.hpp | 3 ++- 3 files changed, 10 insertions(+), 12 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index af66ef6e..89228bd6 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -342,8 +342,11 @@ int main(int argc, char *argv[]) args.AddOption(&hyperreductionSamplingType, "-hrsamptype", "--hrsamplingtype", "Sampling type for the hyperreduction."); args.AddOption(&romOptions.maxNNLSnnz, "-maxnnls", "--max-nnls", - "Maximum nnz for NNLS"); - args.Parse(); + "Maximum number of nonzeros in NNLS solution."); + args.AddOption(&romOptions.tolNNLS, "-tolnnls", "--tol-nnls", + "NNLS solver error tolerance."); + + args.Parse(); if (!args.Good()) { if (mpi.Root()) { diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index ed1438f8..9d94abc8 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -683,11 +683,8 @@ void ROM_Sampler::SetupEQP_Force_Eq(const CAROM::Matrix* snapX, w((i*nqe) + j) = w_el[j]; } - // TODO: input these NNLS parameters? - double tolNNLS = 1.0e-14; - CAROM::Vector sol(ne * nqe, true); - SolveNNLS(rank, tolNNLS, input.maxNNLSnnz, w, Gt, sol); + SolveNNLS(rank, input.tolNNLS, input.maxNNLSnnz, w, Gt, sol); const std::string varName = equationE ? "E" : "V"; WriteSolutionNNLS(sol, "run/nnls" + varName + std::to_string(input.window) + "_" + @@ -971,11 +968,8 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, } } - // TODO: input these NNLS parameters? - double tolNNLS = 1.0e-14; - - CAROM::Vector sol(ne * nqe, true); - SolveNNLS(rank, tolNNLS, input.maxNNLSnnz, w, Gt, sol); + CAROM::Vector sol(ne * nqe, true); + SolveNNLS(rank, input.tolNNLS, input.maxNNLSnnz, w, Gt, sol); const std::string varName = "EC"; // energy conserving WriteSolutionNNLS(sol, "run/nnls" + varName + std::to_string(input.window) diff --git a/rom/laghos_rom.hpp b/rom/laghos_rom.hpp index 87fe723a..9bfa783b 100644 --- a/rom/laghos_rom.hpp +++ b/rom/laghos_rom.hpp @@ -259,7 +259,8 @@ struct ROM_Options bool VTos = false; - int maxNNLSnnz = 0; + int maxNNLSnnz = 0; // max number of NNLS solution nonzeros + double tolNNLS = 1.0e-14; // NNLS solver error tolerance }; static double* getGreedyParam(ROM_Options& romOptions, const char* greedyParam) From fa591bd1e6742952d6f74ed94553a1b556e8550d Mon Sep 17 00:00:00 2001 From: Chris Val Date: Tue, 18 Jul 2023 10:40:11 -0700 Subject: [PATCH 16/39] add gitignore file --- rom/.gitignore | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 rom/.gitignore diff --git a/rom/.gitignore b/rom/.gitignore new file mode 100644 index 00000000..437415e7 --- /dev/null +++ b/rom/.gitignore @@ -0,0 +1,7 @@ +dependencies/ +run/ +laghos +merge +*.o +*.out +job*.sh From 124dd00fff5236b4ed5dd293e17d08d7db80941a Mon Sep 17 00:00:00 2001 From: Chris Val Date: Mon, 31 Jul 2023 16:23:41 -0700 Subject: [PATCH 17/39] use libROM matrix rescaling method --- rom/.gitignore | 1 + rom/laghos_rom.cpp | 30 +++--------------------------- 2 files changed, 4 insertions(+), 27 deletions(-) diff --git a/rom/.gitignore b/rom/.gitignore index 437415e7..ab1ebe02 100644 --- a/rom/.gitignore +++ b/rom/.gitignore @@ -5,3 +5,4 @@ merge *.o *.out job*.sh +logfiles/ diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 9d94abc8..1b122a2e 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -931,8 +931,9 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, } // j -- basis vector } // i -- snapshot - // Rescale every column (NNLS equation) by its max absolute value - //Gt.rescale_cols_max(); + // Rescale every Gt column (NNLS equation) by its max absolute value. + // It seems to help the NNLS solver significantly. + Gt.rescale_cols_max(); CAROM::Vector w(ne * nqe, true); @@ -943,31 +944,6 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX, w((i*nqe) + j) = w_el[j]; } - if (rank == 0) - { - int nGtcols = (NBv + NBe) * nsnap; - int nGtrows = NQ; - for (int jj = 0; jj < nGtcols; jj++) - { - double tmpmax = fabs(Gt(0, jj)); - for (int ii = 1; ii < nGtrows; ii++) - { - if (fabs(Gt(ii, jj)) > tmpmax) - { - tmpmax = fabs(Gt(ii, jj)); - } - } - - if (tmpmax > 1.0e-14) - { - for (int ii = 0; ii < nGtrows; ii++) - { - Gt(ii, jj) /= tmpmax; - } - } - } - } - CAROM::Vector sol(ne * nqe, true); SolveNNLS(rank, input.tolNNLS, input.maxNNLSnnz, w, Gt, sol); From 3fa45efa9d5a181bdfa618b3f9ee6b332598e525 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Tue, 1 Aug 2023 10:04:38 -0700 Subject: [PATCH 18/39] sampling frequency; energy calc every dt --- rom/laghos.cpp | 54 +++++++++++++++++++++++++++++++++++++++++----- rom/laghos_rom.hpp | 3 +++ 2 files changed, 52 insertions(+), 5 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 89228bd6..c174c362 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -345,6 +345,8 @@ int main(int argc, char *argv[]) "Maximum number of nonzeros in NNLS solution."); args.AddOption(&romOptions.tolNNLS, "-tolnnls", "--tol-nnls", "NNLS solver error tolerance."); + args.AddOption(&romOptions.sampfreq, "-sampfreq", "--samp-freq", + "Snapshot sampling frequency."); args.Parse(); if (!args.Good()) @@ -1178,6 +1180,10 @@ int main(int argc, char *argv[]) sampler->SampleSolution(0, 0, (problem == 7) ? 0.0 : -1.0, *S); } samplerTimer.Stop(); + + if (myid == 0) + cout << "Sampling every " << romOptions.sampfreq << + " timestep(s)." << endl; } if (outputTimes) @@ -1820,7 +1826,8 @@ int main(int argc, char *argv[]) } else { - sampler->SampleSolution(t, last_dt, real_pd, *S); + if (unique_steps % romOptions.sampfreq == 0) + sampler->SampleSolution(t, last_dt, real_pd, *S); } bool endWindow = false; @@ -1986,6 +1993,22 @@ int main(int argc, char *argv[]) if (rom_online) { + if (romOptions.hyperreductionSamplingType == eqp_energy) + { + double energy_total, energy_diff; + energy_total = oper->InternalEnergy(*e_gf) + + oper->KineticEnergy(*v_gf); + energy_diff = energy_total - energy_init; + + if (mpi.Root()) + { + cout << "\tE_tot = " << scientific << setprecision(5) + << energy_total + << ",\tE_diff = " << scientific << setprecision(5) + << energy_diff << endl; + } + } + double window_par = t; if (problem == 7) { @@ -2049,10 +2072,31 @@ int main(int argc, char *argv[]) { basis[romOptions.window]->ProjectFOMtoROM(*S, romS); } - if (myid == 0) - { - cout << "Window " << romOptions.window << ": initial romS norm " << romS.Norml2() << endl; - } + if (myid == 0) + { + cout << "Window " << romOptions.window << ": initial romS norm " << romS.Norml2() << endl; + } + + // Recompute and print the energy information after the + // change of basis, before any more timestepping takes + // place. + if (romOptions.hyperreductionSamplingType == eqp_energy) + { + basis[romOptions.window]->LiftROMtoFOM(romS, *S); + + double energy_total, energy_diff; + energy_total = oper->InternalEnergy(*e_gf) + + oper->KineticEnergy(*v_gf); + energy_diff = energy_total - energy_init; + + if (mpi.Root()) + { + cout << "\tE_tot = " << scientific << setprecision(5) + << energy_total + << ",\tE_diff = " << scientific << setprecision(5) + << energy_diff << endl; + } + } delete romOper[romOptions.window-1]; diff --git a/rom/laghos_rom.hpp b/rom/laghos_rom.hpp index 9bfa783b..c7376bb5 100644 --- a/rom/laghos_rom.hpp +++ b/rom/laghos_rom.hpp @@ -261,6 +261,9 @@ struct ROM_Options int maxNNLSnnz = 0; // max number of NNLS solution nonzeros double tolNNLS = 1.0e-14; // NNLS solver error tolerance + + // snapshot sampling frequency (sample every sampfreq timestep) + int sampfreq = 1; }; static double* getGreedyParam(ROM_Options& romOptions, const char* greedyParam) From 4105bae1f856a1496eaccd47ec49ba4f441b0d41 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Thu, 3 Aug 2023 17:15:59 -0700 Subject: [PATCH 19/39] restore sampling call flags --- rom/laghos_rom.cpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 1b122a2e..10d044cf 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -57,8 +57,7 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p SetStateVariables(S); SetStateVariableRates(dt); - const bool sampleX = (hyperreductionSamplingType == eqp_energy) ? true : - generator_X->isNextSample(t); + const bool sampleX = generator_X->isNextSample(t); Vector dSdt; if (!sns && rhsBasis) @@ -99,8 +98,7 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p } } - const bool sampleV = (hyperreductionSamplingType == eqp_energy) ? true : - generator_V->isNextSample(t); + const bool sampleV = generator_V->isNextSample(t); //TODO: use this, plus generator_Fv->computeNextSampleTime? So far, it seems sampleV == true on every step. //const bool sampleFv = generator_Fv->isNextSample(t); @@ -154,8 +152,7 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p } } - const bool sampleE = (hyperreductionSamplingType == eqp_energy) ? true : - generator_E->isNextSample(t); + const bool sampleE = generator_E->isNextSample(t); if (sampleE) { From 8f5ddc398d844d8cbda1df641b54498da2a6570b Mon Sep 17 00:00:00 2001 From: Chris Val Date: Fri, 4 Aug 2023 10:50:54 -0700 Subject: [PATCH 20/39] minor matrix allocation change --- rom/laghos_rom.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 10d044cf..869350d7 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -1104,7 +1104,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) // increased by 1 to accomodate the energy identity. CAROM::Matrix *tBasisV = basisV->getFirstNColumns(cutoff[1]); - CAROM::Matrix *tBasisE = new CAROM::Matrix(tL2size, cutoff[2], false); + CAROM::Matrix *tBasisE = new CAROM::Matrix(tL2size, cutoff[2], true); // Get the first cutoff[2]-1 columns of basisE for (int i = 0; i < tL2size; i++) From cd57d9d6daf5d50b011d76715060b4533d2d93e5 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Fri, 11 Aug 2023 13:07:32 -0700 Subject: [PATCH 21/39] fix online loop timing --- rom/laghos.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index c174c362..c95624a7 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -1993,8 +1993,9 @@ int main(int argc, char *argv[]) if (rom_online) { - if (romOptions.hyperreductionSamplingType == eqp_energy) + if (romOptions.hyperreductionSamplingType == eqp_energy) { + timeLoopTimer.Stop(); double energy_total, energy_diff; energy_total = oper->InternalEnergy(*e_gf) + oper->KineticEnergy(*v_gf); @@ -2007,6 +2008,7 @@ int main(int argc, char *argv[]) << ",\tE_diff = " << scientific << setprecision(5) << energy_diff << endl; } + timeLoopTimer.Start(); } double window_par = t; @@ -2082,6 +2084,7 @@ int main(int argc, char *argv[]) // place. if (romOptions.hyperreductionSamplingType == eqp_energy) { + timeLoopTimer.Stop(); basis[romOptions.window]->LiftROMtoFOM(romS, *S); double energy_total, energy_diff; @@ -2096,6 +2099,7 @@ int main(int argc, char *argv[]) << ",\tE_diff = " << scientific << setprecision(5) << energy_diff << endl; } + timeLoopTimer.Start(); } delete romOper[romOptions.window-1]; @@ -2158,6 +2162,7 @@ int main(int argc, char *argv[]) if (romOptions.hyperreductionSamplingType == eqp_energy) { + timeLoopTimer.Stop(); double energy_total, energy_diff; energy_total = oper->InternalEnergy(*e_gf) + oper->KineticEnergy(*v_gf); @@ -2170,6 +2175,7 @@ int main(int argc, char *argv[]) << ",\tE_diff = " << scientific << setprecision(5) << energy_diff << endl; } + timeLoopTimer.Start(); } // Make sure all ranks have sent their 'v' solution before initiating From eda3d08690625f90ae66a339420975d863f94865 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Mon, 4 Sep 2023 11:19:48 -0400 Subject: [PATCH 22/39] windowing: online stage changes --- rom/laghos.cpp | 8 ++ rom/laghos_rom.cpp | 188 +++++++++++++++++++++++++++++++++++++++------ rom/laghos_rom.hpp | 3 + 3 files changed, 175 insertions(+), 24 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index c95624a7..67c731b2 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -2072,6 +2072,14 @@ int main(int argc, char *argv[]) if (!romOptions.use_sample_mesh) { + if (romOptions.hyperreductionSamplingType == eqp_energy) + { + // Add the corresponding lifted solution vector + // to the velocity and energy bases to ensure + // energy is conserved across windows. + basis[romOptions.window]->AddLastCol_V(*S); + basis[romOptions.window]->AddLastCol_E(*S); + } basis[romOptions.window]->ProjectFOMtoROM(*S, romS); } if (myid == 0) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 869350d7..de3e0b56 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -1275,21 +1275,43 @@ ROM_Basis::ROM_Basis(ROM_Options const& input, MPI_Comm comm_, const double sFac mfH1.SetSize(tH1size); mfL2.SetSize(tL2size); - if (hyperreductionSamplingType == eqp_energy) - basisE = new CAROM::Matrix(tL2size, rdime, true); - - ReadSolutionBases(input.window); - + // This code block is used to set the sizes of the basis matrices + // correctly, before the reading of the actual data takes place + // within ReadSolutionBases following right after. + // In this way, once the basis data has been read, the matrices + // can be extended by adding the required column vectors. if (hyperreductionSamplingType == eqp_energy) { - // Form the energy identity and include it as the last basis vector. - Vector unitE(tL2size); - unitE = 1.0; - for (int i = 0; i < tL2size; i++) - (*basisE)(i, rdime-1) = unitE[i]; - basisE->orthogonalize(); + if (input.window == 0) + { + // For the first window, only the energy basis need be extended. + + // The energy basis is extended by adding the energy identity. + // Dimension rdime has already been increased by 1 in the + // offline stage + basisE = new CAROM::Matrix(tL2size, rdime, true); + } + else + { + // For any window other than the first, both energy and + // velocity bases need be extended. + + // The velocity basis is extended by adding the current lifted + // velocity solution vector. + // Dimension rdimv has already been increased by 1 in the + // offline stage. + basisV = new CAROM::Matrix(tH1size, rdimv, true); + + // The energy basis is extended by adding the energy identity + // and the current lifted energy solution vector. + // Dimension rdime has already been increased by 2 in the + // offline stage. + basisE = new CAROM::Matrix(tL2size, rdime, true); + } } + ReadSolutionBases(input.window); + if (spaceTime) { ReadTemporalBases(input.window); @@ -2207,24 +2229,110 @@ int ROM_Basis::SolutionSizeFOM() const void ROM_Basis::ReadSolutionBases(const int window) { - if (!useVX) - basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + std::to_string(window) + basisIdentifier, tH1size, rdimv); - - // In the energy-conserving EQP case we read the first rdime-1 basis - // vectors, since the rdime parameter has been increased by 1 to - // accommodate the addition of the energy identity. + // Velocity basis formation. if (hyperreductionSamplingType == eqp_energy) { - int tmp_rdime = rdime - 1; - CAROM::Matrix *tmp_basisE = 0; + if (window == 0) + { + // For the first window, read all rdimv basis vectors normally. + basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + + std::to_string(window) + basisIdentifier, tH1size, rdimv); + } + else + { + // For any window other than the first, read the first rdimv-1 + // basis vectors, since rdimv has been increased by 1 to accomodate + // the addition of the lifted velocity solution vector. + int tmp_rdimv = rdimv - 1; + CAROM::Matrix *tmp_basisV = 0; + + tmp_basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + + std::to_string(window) + basisIdentifier, tH1size, tmp_rdimv); + + for (int i = 0; i < tH1size; i++) + { + for (int j = 0; j < tmp_rdimv; j++) + (*basisV)(i, j) = (*tmp_basisV)(i, j); + } + delete tmp_basisV; + + // The addition of the lifted velocity solution vector as the last + // column vector takes place later in the main driver "laghos.cpp" + // at the time of the window change, because the current solution + // vector S must be known. + } + } + else + { + if (!useVX) + { + basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + + std::to_string(window) + basisIdentifier, tH1size, rdimv); + } + } - tmp_basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + - std::to_string(window) + basisIdentifier, tL2size, tmp_rdime); + // Energy basis formation. + if (hyperreductionSamplingType == eqp_energy) + { + if (window == 0) + { + // For the first window, read the first rdime-1 basis vectors, + // since rdime has been increased by 1 to accomodate the addition + // of the energy identity. + int tmp_rdime = rdime - 1; + CAROM::Matrix *tmp_basisE = 0; + + tmp_basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + + std::to_string(window) + basisIdentifier, tL2size, tmp_rdime); + + for (int i = 0; i < tL2size; i++) + { + for (int j = 0; j < tmp_rdime; j++) + (*basisE)(i, j) = (*tmp_basisE)(i, j); + } + delete tmp_basisE; - for (int i = 0; i < tL2size; i++) + // Add the energy identity as the last column vector. + Vector unitE(tL2size); + unitE = 1.0; + for (int i = 0; i < tL2size; i++) + { + (*basisE)(i, tmp_rdime) = unitE[i]; + } + basisE->orthogonalize(); + } + else { - for (int j = 0; j < tmp_rdime; j++) - (*basisE)(i, j) = (*tmp_basisE)(i, j); + // For any window other than the first, read the first rdime-2 + // basis vectors, since rdime has been increased by 2 to accomodate + // the addition of the energy identity and the lifted energy + // solution vector. + int tmp_rdime = rdime - 2; + CAROM::Matrix *tmp_basisE = 0; + + tmp_basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + + std::to_string(window) + basisIdentifier, tL2size, tmp_rdime); + + for (int i = 0; i < tL2size; i++) + { + for (int j = 0; j < tmp_rdime; j++) + (*basisE)(i, j) = (*tmp_basisE)(i, j); + } + delete tmp_basisE; + + // Add the energy identity as the penultimate column vector. + Vector unitE(tL2size); + unitE = 1.0; + for (int i = 0; i < tL2size; i++) + { + (*basisE)(i, tmp_rdime) = unitE[i]; + } + basisE->orthogonalize(); + + // The addition of the lifted energy solution vector as the last + // column vector takes place later in the main driver "laghos.cpp" + // at the time of the window change, because the current solution + // vector S must be known. } } else @@ -2404,6 +2512,38 @@ void ROM_Basis::ProjectFOMtoROM_V(Vector const& f, Vector & r, const bool timeDe r[i] = (*rV)(i); } +// f is a full vector, not a true vector +void ROM_Basis::AddLastCol_V(Vector const& f) +{ + MFEM_VERIFY(f.Size() == (2*H1size) + L2size, ""); + + for (int i=0; iGetTrueDofs(mfH1); + + for (int i = 0; i < tH1size; i++) + (*basisV)(i, rdimv-1) = mfH1[i]; + + basisV->orthogonalize(); +} + +// f is a full vector, not a true vector +void ROM_Basis::AddLastCol_E(Vector const& f) +{ + MFEM_VERIFY(f.Size() == (2*H1size) + L2size, ""); + + for (int i=0; iGetTrueDofs(mfL2); + + for (int i = 0; i < tL2size; i++) + (*basisE)(i, rdime-1) = mfL2[i]; + + basisE->orthogonalize(); +} + // f is a full vector, not a true vector void ROM_Basis::LiftROMtoFOM(Vector const& r, Vector & f) { diff --git a/rom/laghos_rom.hpp b/rom/laghos_rom.hpp index c7376bb5..634e1222 100644 --- a/rom/laghos_rom.hpp +++ b/rom/laghos_rom.hpp @@ -858,6 +858,9 @@ class ROM_Basis void ProjectFOMtoROM_V(Vector const& f, Vector & r, const bool timeDerivative=false); + void AddLastCol_V(Vector const& f); + void AddLastCol_E(Vector const& f); + void LiftROMtoFOM(Vector const& r, Vector & f); void LiftROMtoFOM_dVdt(Vector const& r, Vector & f); void LiftROMtoFOM_dEdt(Vector const& r, Vector & f); From 5e95bae104ea649aa694b64cfc16d40fe5ee409d Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Mon, 4 Sep 2023 11:51:36 -0400 Subject: [PATCH 23/39] windowing: offline stage changes --- rom/laghos_rom.cpp | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index de3e0b56..44b0c1ef 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -1087,10 +1087,8 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) { if (rank == 0) { - // For the energy-conserving EQP case, increase the energy basis - // dimension by 1 to accomodate the addition of the energy - // identity. - + // Increase the energy basis dimension by 1 to accomodate the + // addition of the energy identity. cutoff[2] += 1; } @@ -1103,7 +1101,6 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) // For the energy basis, the cutoff parameter has already been // increased by 1 to accomodate the energy identity. CAROM::Matrix *tBasisV = basisV->getFirstNColumns(cutoff[1]); - CAROM::Matrix *tBasisE = new CAROM::Matrix(tL2size, cutoff[2], true); // Get the first cutoff[2]-1 columns of basisE @@ -1116,8 +1113,10 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) // Form the energy identity and include it as the last basis vector. Vector unitE(tL2size); unitE = 1.0; + for (int i = 0; i < tL2size; i++) (*tBasisE)(i, cutoff[2]-1) = unitE[i]; + tBasisE->orthogonalize(); SetupEQP_Force(generator_X->getSnapshotMatrix(), @@ -1127,6 +1126,17 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) delete tBasisV; delete tBasisE; + + if (rank == 0 && input.window > 0) + { + // If not first window, increase the dimension of the velocity and + // energy bases by 1, to accomodate the future addition of the + // lifted solution vector in the online stage. + // The change is made here so that the right dimensions are written + // in the window parameters file in the caller's scope. + cutoff[1] += 1; + cutoff[2] += 1; + } } delete generator_X; From 85948c9a7f2c501ad40861f721df78175e88e0fc Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Tue, 26 Sep 2023 09:40:58 -0700 Subject: [PATCH 24/39] form ROM basis/operator at time of window change --- rom/laghos.cpp | 78 ++++++++++++++++++++++++++++++++++---------------- 1 file changed, 53 insertions(+), 25 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 67c731b2..2a167479 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -1244,18 +1244,37 @@ int main(int argc, char *argv[]) if (dtc > 0.0) dt = dtc; if (usingWindows) { - // Construct the ROM_Basis for each window. - for (romOptions.window = numWindows-1; romOptions.window >= 0; --romOptions.window) + if (romOptions.hyperreductionSamplingType != eqp_energy) { + // Construct the ROM_Basis for each window. + for (romOptions.window = numWindows-1; romOptions.window >= 0; --romOptions.window) + { + SetWindowParameters(twparam, romOptions); + basis[romOptions.window] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV); + if (!romOptions.hyperreduce_prep) + { + romOper[romOptions.window] = new ROM_Operator(romOptions, basis[romOptions.window], rho_coeff, mat_coeff, order_e, source, + visc, vort, cfl, p_assembly, cg_tol, cg_max_iter, ftz_tol, &H1FEC, &L2FEC); + } + } + romOptions.window = 0; + } + else + { + // Construct the ROM basis and operator only for the first + // window; the ones for the next windows will be constructed + // at the time of each window change, because the current + // vector S is needed for their construction. + // In particular, forming the bases and reduced inverse mass + // matrices requires knoweldege of vector S. SetWindowParameters(twparam, romOptions); - basis[romOptions.window] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV); + basis[0] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV, ×teps); if (!romOptions.hyperreduce_prep) { - romOper[romOptions.window] = new ROM_Operator(romOptions, basis[romOptions.window], rho_coeff, mat_coeff, order_e, source, - visc, vort, cfl, p_assembly, cg_tol, cg_max_iter, ftz_tol, &H1FEC, &L2FEC); + romOper[0] = new ROM_Operator(romOptions, basis[0], rho_coeff, mat_coeff, order_e, source, visc, vort, cfl, p_assembly, + cg_tol, cg_max_iter, ftz_tol, &H1FEC, &L2FEC, ×teps); } } - romOptions.window = 0; } else { @@ -1590,12 +1609,6 @@ int main(int argc, char *argv[]) romOper[0]->ApplyHyperreduction(romS); } - // // TODO: do we want that for the energy-conserving EQP? - // if (rom_online && romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) - // { - // romOper[0]->ApplyHyperreduction(romS); - // } - double windowEndpoint = 0.0; double windowOverlapMidpoint = 0.0; for (ti = 1; !last_step; ti++) @@ -2045,6 +2058,34 @@ int main(int argc, char *argv[]) int rdimeprev = romOptions.dimE; SetWindowParameters(twparam, romOptions); + if (romOptions.hyperreductionSamplingType == eqp_energy) + { + timeLoopTimer.Stop(); + onlinePreprocessTimer.Start(); + + // Form the ROM basis and operator for the new window. + basis[romOptions.window] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV, ×teps); + + onlinePreprocessTimer.Stop(); + timeLoopTimer.Start(); + + // Add the last lifted solution vector as the last + // column in the bases. + basis[romOptions.window]->AddLastCol_V(*S); + basis[romOptions.window]->AddLastCol_E(*S); + + timeLoopTimer.Stop(); + onlinePreprocessTimer.Start(); + + if (!romOptions.hyperreduce_prep) + { + romOper[romOptions.window] = new ROM_Operator(romOptions, basis[romOptions.window], rho_coeff, mat_coeff, + order_e, source, visc, vort, cfl, p_assembly, + cg_tol, cg_max_iter, ftz_tol, &H1FEC, &L2FEC, ×teps); + } + onlinePreprocessTimer.Stop(); + timeLoopTimer.Start(); + } if (romOptions.use_sample_mesh) { basis[romOptions.window]->ProjectFromPreviousWindow(romOptions, romS, romOptions.window, rdimxprev, rdimvprev, rdimeprev); @@ -2072,14 +2113,6 @@ int main(int argc, char *argv[]) if (!romOptions.use_sample_mesh) { - if (romOptions.hyperreductionSamplingType == eqp_energy) - { - // Add the corresponding lifted solution vector - // to the velocity and energy bases to ensure - // energy is conserved across windows. - basis[romOptions.window]->AddLastCol_V(*S); - basis[romOptions.window]->AddLastCol_E(*S); - } basis[romOptions.window]->ProjectFOMtoROM(*S, romS); } if (myid == 0) @@ -2117,11 +2150,6 @@ int main(int argc, char *argv[]) romOper[romOptions.window]->ApplyHyperreduction(romS); } - //if (romOptions.hyperreduce && romOptions.hyperreductionSamplingType == eqp_energy) - //{ - // romOper[romOptions.window]->ApplyHyperreduction(romS); - //} - if (problem == 7 && romOptions.indicatorType == penetrationDistance) { std::string pd_weight_outputPath = testing_parameter_outputPath + "/pd_weight" + to_string(romOptions.window); From 593497e8f3a4e4a84ce12bc22fc37125793819a1 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Tue, 26 Sep 2023 09:42:44 -0700 Subject: [PATCH 25/39] fix basis formation in restore phase --- rom/laghos.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 2a167479..5a49d0c6 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -1536,7 +1536,16 @@ int main(int argc, char *argv[]) SetWindowParameters(twparam, romOptions); basis[romOptions.window-1]->LiftROMtoFOM(romS, *S); delete basis[romOptions.window-1]; + basis[romOptions.window] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV); + + if (romOptions.hyperreductionSamplingType == eqp_energy) + { + // Add the last lifted solution vector as the last + // column in the bases. + basis[romOptions.window]->AddLastCol_V(*S); + basis[romOptions.window]->AddLastCol_E(*S); + } basis[romOptions.window]->Init(romOptions, *S); if (romOptions.mergeXV) From 09e72685544aa2ac89033e8485787c5dc54850bc Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Fri, 20 Oct 2023 09:14:48 -0700 Subject: [PATCH 26/39] add Dylan's EQP bug fixes --- rom/laghos.cpp | 8 ++++---- rom/laghos_rom.cpp | 13 +++++++++---- rom/laghos_rom.hpp | 7 ++++--- 3 files changed, 17 insertions(+), 11 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 5a49d0c6..e9ec5c08 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -1914,7 +1914,7 @@ int main(int argc, char *argv[]) if (samplerLast->MaxNumSamples() >= windowNumSamples + windowOverlapSamples || last_step) { - samplerLast->Finalize(cutoff, romOptions); + samplerLast->Finalize(cutoff, romOptions, *S); if (last_step) { // Let samplerLast define the final window, discarding the sampler window. @@ -1960,7 +1960,7 @@ int main(int argc, char *argv[]) else { // Form the reduced bases from the snapshots. - sampler->Finalize(cutoff, romOptions); + sampler->Finalize(cutoff, romOptions, *S); } if (myid == 0 && romOptions.parameterID == -1) { outfile_twp << windowEndpoint << ", " << cutoff[0] << ", " << cutoff[1] << ", " << cutoff[2]; @@ -2323,9 +2323,9 @@ int main(int argc, char *argv[]) else { if (samplerLast) - samplerLast->Finalize(cutoff, romOptions); + samplerLast->Finalize(cutoff, romOptions, *S); else if (sampler) - sampler->Finalize(cutoff, romOptions); + sampler->Finalize(cutoff, romOptions, *S); } basisConstructionTimer.Stop(); diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 44b0c1ef..039539c7 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -954,7 +954,8 @@ void ROM_Sampler::SetupEQP_Force(const CAROM::Matrix* snapX, const CAROM::Matrix* snapE, const CAROM::Matrix* basisV, const CAROM::Matrix* basisE, - ROM_Options const& input) + ROM_Options const& input, + Vector const& sol) { MFEM_VERIFY(basisV->numRows() == input.H1FESpace->GetTrueVSize(), ""); MFEM_VERIFY(basisE->numRows() == input.L2FESpace->GetTrueVSize(), ""); @@ -974,9 +975,13 @@ void ROM_Sampler::SetupEQP_Force(const CAROM::Matrix* snapX, // energy-conserving EQP: one combined rule for velocity and energy SetupEQP_En_Force_Eq(snapX, snapV, snapE, basisV, basisE, input); } + + // Call this to call UpdateQuadratureData and restore the FOM state. + input.FOMoper->GetTimeStepEstimate(sol); } -void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) +void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, + Vector const& sol) { if (writeSnapshots) { @@ -1077,7 +1082,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) SetupEQP_Force(generator_X->getSnapshotMatrix(), generator_V->getSnapshotMatrix(), generator_E->getSnapshotMatrix(), - tBasisV, tBasisE, input); + tBasisV, tBasisE, input, sol); delete tBasisV; delete tBasisE; @@ -1122,7 +1127,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input) SetupEQP_Force(generator_X->getSnapshotMatrix(), generator_V->getSnapshotMatrix(), generator_E->getSnapshotMatrix(), - tBasisV, tBasisE, input); + tBasisV, tBasisE, input, sol); delete tBasisV; delete tBasisE; diff --git a/rom/laghos_rom.hpp b/rom/laghos_rom.hpp index 634e1222..9c9b1c7c 100644 --- a/rom/laghos_rom.hpp +++ b/rom/laghos_rom.hpp @@ -492,7 +492,7 @@ class ROM_Sampler // TODO: update the following comment, since there should now be a maximum of 1 time interval now. const int max_model_dim_est = int(input.t_final/input.initial_dt + 0.5) + 100; // Note that this is a rough estimate which may be exceeded, resulting in multiple libROM basis time intervals. - const int max_model_dim = (input.max_dim > 0) ? input.max_dim : max_model_dim_est; + const int max_model_dim = (input.max_dim > 0) ? input.max_dim + 1 : max_model_dim_est; std::cout << rank << ": max_model_dim " << max_model_dim << std::endl; @@ -632,7 +632,7 @@ class ROM_Sampler void SampleSolution(const double t, const double dt, const double pd, Vector const& S); - void Finalize(Array &cutoff, ROM_Options& input); + void Finalize(Array &cutoff, ROM_Options& input, Vector const& sol); int MaxNumSamples() { @@ -826,7 +826,8 @@ class ROM_Sampler } void SetupEQP_Force(const CAROM::Matrix* snapX, const CAROM::Matrix* snapV, const CAROM::Matrix* snapE, - const CAROM::Matrix* basisV, const CAROM::Matrix* basisE, ROM_Options const& input); + const CAROM::Matrix* basisV, const CAROM::Matrix* basisE, ROM_Options const& input, + Vector const& sol); void SetupEQP_Force_Eq(const CAROM::Matrix* snapX, const CAROM::Matrix* snapV, const CAROM::Matrix* snapE, const CAROM::Matrix* basisV, const CAROM::Matrix* basisE, ROM_Options const& input, From b0e6d08428b1f7a352ac0b8950e84d3c00431317 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Fri, 3 Nov 2023 11:33:21 -0700 Subject: [PATCH 27/39] fix L2 errors calculation --- rom/laghos.cpp | 57 +++++++++++++++++++++++++++++++------------------- 1 file changed, 35 insertions(+), 22 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index e9ec5c08..96dfa2f7 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -2129,28 +2129,28 @@ int main(int argc, char *argv[]) cout << "Window " << romOptions.window << ": initial romS norm " << romS.Norml2() << endl; } - // Recompute and print the energy information after the - // change of basis, before any more timestepping takes - // place. - if (romOptions.hyperreductionSamplingType == eqp_energy) - { - timeLoopTimer.Stop(); - basis[romOptions.window]->LiftROMtoFOM(romS, *S); - - double energy_total, energy_diff; - energy_total = oper->InternalEnergy(*e_gf) + - oper->KineticEnergy(*v_gf); - energy_diff = energy_total - energy_init; - - if (mpi.Root()) - { - cout << "\tE_tot = " << scientific << setprecision(5) - << energy_total - << ",\tE_diff = " << scientific << setprecision(5) - << energy_diff << endl; - } - timeLoopTimer.Start(); - } + //if (romOptions.hyperreductionSamplingType == eqp_energy) + //{ + // // Recompute and print the energy information after the + // // change of basis, before any more timestepping takes + // // place. + // timeLoopTimer.Stop(); + // basis[romOptions.window]->LiftROMtoFOM(romS, *S); + + // double energy_total, energy_diff; + // energy_total = oper->InternalEnergy(*e_gf) + + // oper->KineticEnergy(*v_gf); + // energy_diff = energy_total - energy_init; + + // if (mpi.Root()) + // { + // cout << "\tE_tot = " << scientific << setprecision(5) + // << energy_total + // << ",\tE_diff = " << scientific << setprecision(5) + // << energy_diff << endl; + // } + // timeLoopTimer.Start(); + //} delete romOper[romOptions.window-1]; @@ -2297,6 +2297,19 @@ int main(int argc, char *argv[]) outfile_tw_steps.close(); } + if (rom_online && romOptions.hyperreductionSamplingType == eqp_energy) + { + if (myid == 0) cout << endl; + if (myid == 0) cout << "Solution errors calculation." << endl; + + PrintDiffParGridFunction(normtype, myid, testing_parameter_outputPath + + "/Sol_Position" + romOptions.basisIdentifier, x_gf); + PrintDiffParGridFunction(normtype, myid, testing_parameter_outputPath + + "/Sol_Velocity" + romOptions.basisIdentifier, v_gf); + PrintDiffParGridFunction(normtype, myid, testing_parameter_outputPath + + "/Sol_Energy" + romOptions.basisIdentifier, e_gf); + } + if (romOptions.use_sample_mesh) { if (romOptions.GramSchmidt && !spaceTime) From 473c0e3af28784fd6a7c6035a25ff02ac7e22ef5 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Fri, 3 Nov 2023 11:35:27 -0700 Subject: [PATCH 28/39] use incremental GS method --- rom/laghos_rom.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 039539c7..23f47a38 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -1122,7 +1122,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, for (int i = 0; i < tL2size; i++) (*tBasisE)(i, cutoff[2]-1) = unitE[i]; - tBasisE->orthogonalize(); + tBasisE->orthogonalize_last(); SetupEQP_Force(generator_X->getSnapshotMatrix(), generator_V->getSnapshotMatrix(), @@ -2314,7 +2314,7 @@ void ROM_Basis::ReadSolutionBases(const int window) { (*basisE)(i, tmp_rdime) = unitE[i]; } - basisE->orthogonalize(); + basisE->orthogonalize_last(); } else { @@ -2342,7 +2342,7 @@ void ROM_Basis::ReadSolutionBases(const int window) { (*basisE)(i, tmp_rdime) = unitE[i]; } - basisE->orthogonalize(); + basisE->orthogonalize_last(rdime-1); // The addition of the lifted energy solution vector as the last // column vector takes place later in the main driver "laghos.cpp" @@ -2540,7 +2540,7 @@ void ROM_Basis::AddLastCol_V(Vector const& f) for (int i = 0; i < tH1size; i++) (*basisV)(i, rdimv-1) = mfH1[i]; - basisV->orthogonalize(); + basisV->orthogonalize_last(); } // f is a full vector, not a true vector @@ -2556,7 +2556,7 @@ void ROM_Basis::AddLastCol_E(Vector const& f) for (int i = 0; i < tL2size; i++) (*basisE)(i, rdime-1) = mfL2[i]; - basisE->orthogonalize(); + basisE->orthogonalize_last(); } // f is a full vector, not a true vector From 242155c765def5341d4284d769e608144f4271e1 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Fri, 3 Nov 2023 13:38:27 -0700 Subject: [PATCH 29/39] use snapshots in EQP derivation --- rom/laghos_rom.cpp | 134 +++++++++++++++++++++++++++++++-------------- 1 file changed, 92 insertions(+), 42 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 23f47a38..0d254baf 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -1090,58 +1090,108 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, if (input.hyperreductionSamplingType == eqp_energy) { - if (rank == 0) + if (rank == 0 && input.window == 0) { - // Increase the energy basis dimension by 1 to accomodate the - // addition of the energy identity. + // If first window, increase the energy basis dimension by 1 + // to accomodate the energy identity. + // No change in the velocity bases dimension. cutoff[2] += 1; } + else if (rank == 0 && input.window > 0) + { + // If not first window, increase the dimension of the velocity + // basis by 1 to accomodate the velocity solution snapshot, + // and the energy basis dimension by 2 to accomodate the energy + // identity and energy solution snapshot. + // We add the FOM snapshots in the bases so that they are used + // in deriving the EQP rule. + // Also, the change is made here so that the right basis + // dimensions are written in the window parameters file in the + // caller's scope. + cutoff[1] += 1; + cutoff[2] += 2; + } + + MPI_Bcast(cutoff.GetData(), cutoff.Size(), MPI_INT, 0, MPI_COMM_WORLD); + + const CAROM::Matrix *snapX = generator_X->getSnapshotMatrix(); + const CAROM::Matrix *snapV = generator_V->getSnapshotMatrix(); + const CAROM::Matrix *snapE = generator_E->getSnapshotMatrix(); const CAROM::Matrix *basisV = generator_V->getSpatialBasis(); const CAROM::Matrix *basisE = generator_E->getSpatialBasis(); - MPI_Bcast(cutoff.GetData(), cutoff.Size(), MPI_INT, 0, MPI_COMM_WORLD); - - // Truncate the bases. - // For the energy basis, the cutoff parameter has already been - // increased by 1 to accomodate the energy identity. - CAROM::Matrix *tBasisV = basisV->getFirstNColumns(cutoff[1]); + // Form the reduced bases. + CAROM::Matrix *tBasisV = new CAROM::Matrix(tH1size, cutoff[1], true); CAROM::Matrix *tBasisE = new CAROM::Matrix(tL2size, cutoff[2], true); + + if (input.window == 0) + { + // Get the first cutoff[1] columns of basisV + for (int i = 0; i < tH1size; i++) + for (int j = 0; j < cutoff[1]; j++) + (*tBasisV)(i, j) = (*basisV)(i, j); + + // Get the first cutoff[2] - 1 columns of basisE + for (int i = 0; i < tL2size; i++) + for (int j = 0; j < cutoff[2] - 1; j++) + (*tBasisE)(i, j) = (*basisE)(i, j); + + // Add the energy identity as the last E basis column + // and reorthonormalize. + Vector unitE(tL2size); + unitE = 1.0; + for (int i = 0; i < tL2size; i++) + (*tBasisE)(i, cutoff[2] - 1) = unitE[i]; + + tBasisE->orthogonalize_last(); + } + else if (input.window > 0) + { + // Get the first cutoff[1] - 1 columns of basisV + for (int i = 0; i < tH1size; i++) + for (int j = 0; j < cutoff[1] - 1; j++) + (*tBasisV)(i, j) = (*basisV)(i, j); + + // Add the first V snapshot as the last V basis column + // and reorthonormalize. + // The current window's first snapshot is the same as the + // previous window's last snapshot, since we are not using + // offset vectors. + for (int i = 0; i < tH1size; i++) + (*tBasisV)(i, cutoff[1] - 1) = (*snapV)(i, 0); + + tBasisV->orthogonalize_last(); + + // Get the first cutoff[2] - 2 columns of basisE + for (int i = 0; i < tL2size; i++) + for (int j = 0; j < cutoff[2] - 2; j++) + (*tBasisE)(i, j) = (*basisE)(i, j); + + // Add the energy identity as the penultimate E basis column + // and reorthonormalize. + Vector unitE(tL2size); + unitE = 1.0; + for (int i = 0; i < tL2size; i++) + (*tBasisE)(i, cutoff[2] - 2) = unitE[i]; + + tBasisE->orthogonalize_last(cutoff[2] - 1); + + // Add the first E snapshot as the last E basis column + // and reorthonormalize. + // The current window's first snapshot is the same as the + // previous window's last snapshot, since we are not using + // offset vectors. + for (int i = 0; i < tL2size; i++) + (*tBasisE)(i, cutoff[2] - 1) = (*snapE)(i, 0); + + tBasisE->orthogonalize_last(); + } - // Get the first cutoff[2]-1 columns of basisE - for (int i = 0; i < tL2size; i++) - { - for (int j = 0; j < cutoff[2] - 1; j++) - (*tBasisE)(i, j) = (*basisE)(i, j); - } - - // Form the energy identity and include it as the last basis vector. - Vector unitE(tL2size); - unitE = 1.0; - - for (int i = 0; i < tL2size; i++) - (*tBasisE)(i, cutoff[2]-1) = unitE[i]; - - tBasisE->orthogonalize_last(); - - SetupEQP_Force(generator_X->getSnapshotMatrix(), - generator_V->getSnapshotMatrix(), - generator_E->getSnapshotMatrix(), - tBasisV, tBasisE, input, sol); - - delete tBasisV; - delete tBasisE; + SetupEQP_Force(snapX, snapV, snapE, tBasisV, tBasisE, input, sol); - if (rank == 0 && input.window > 0) - { - // If not first window, increase the dimension of the velocity and - // energy bases by 1, to accomodate the future addition of the - // lifted solution vector in the online stage. - // The change is made here so that the right dimensions are written - // in the window parameters file in the caller's scope. - cutoff[1] += 1; - cutoff[2] += 1; - } + // delete snapX, snapV, snapE; + delete tBasisV, tBasisE; } delete generator_X; From cf55f61c2271a20b3d51eb6af7f12845d63cbfe1 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Thu, 7 Dec 2023 08:08:38 -0800 Subject: [PATCH 30/39] snapshots matrix bug fix --- rom/laghos.cpp | 44 ++++++++++++++--------------- rom/laghos_rom.cpp | 70 +++++++++++++++++++++++----------------------- 2 files changed, 57 insertions(+), 57 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 96dfa2f7..4afabde9 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -2129,28 +2129,28 @@ int main(int argc, char *argv[]) cout << "Window " << romOptions.window << ": initial romS norm " << romS.Norml2() << endl; } - //if (romOptions.hyperreductionSamplingType == eqp_energy) - //{ - // // Recompute and print the energy information after the - // // change of basis, before any more timestepping takes - // // place. - // timeLoopTimer.Stop(); - // basis[romOptions.window]->LiftROMtoFOM(romS, *S); - - // double energy_total, energy_diff; - // energy_total = oper->InternalEnergy(*e_gf) + - // oper->KineticEnergy(*v_gf); - // energy_diff = energy_total - energy_init; - - // if (mpi.Root()) - // { - // cout << "\tE_tot = " << scientific << setprecision(5) - // << energy_total - // << ",\tE_diff = " << scientific << setprecision(5) - // << energy_diff << endl; - // } - // timeLoopTimer.Start(); - //} + if (romOptions.hyperreductionSamplingType == eqp_energy) + { + // Recompute and print the energy information after the + // change of basis, before any more timestepping takes + // place. + timeLoopTimer.Stop(); + basis[romOptions.window]->LiftROMtoFOM(romS, *S); + + double energy_total, energy_diff; + energy_total = oper->InternalEnergy(*e_gf) + + oper->KineticEnergy(*v_gf); + energy_diff = energy_total - energy_init; + + if (mpi.Root()) + { + cout << "\tE_tot = " << scientific << setprecision(5) + << energy_total + << ",\tE_diff = " << scientific << setprecision(5) + << energy_diff << endl; + } + timeLoopTimer.Start(); + } delete romOper[romOptions.window-1]; diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 0d254baf..8dafd4ee 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -983,6 +983,10 @@ void ROM_Sampler::SetupEQP_Force(const CAROM::Matrix* snapX, void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, Vector const& sol) { + CAROM::Matrix Xsnap0(*(generator_X->getSnapshotMatrix())); + CAROM::Matrix Vsnap0(*(generator_V->getSnapshotMatrix())); + CAROM::Matrix Esnap0(*(generator_E->getSnapshotMatrix())); + if (writeSnapshots) { if (!useXV) generator_X->writeSnapshot(); @@ -1088,43 +1092,39 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, delete tBasisE; } - if (input.hyperreductionSamplingType == eqp_energy) - { - if (rank == 0 && input.window == 0) - { + if (input.hyperreductionSamplingType == eqp_energy) + { + if (rank == 0 && input.window == 0) + { // If first window, increase the energy basis dimension by 1 // to accomodate the energy identity. // No change in the velocity bases dimension. - cutoff[2] += 1; - } + cutoff[2] += 1; + } else if (rank == 0 && input.window > 0) - { - // If not first window, increase the dimension of the velocity - // basis by 1 to accomodate the velocity solution snapshot, + { + // If not first window, increase the dimension of the velocity + // basis by 1 to accomodate the velocity solution snapshot, // and the energy basis dimension by 2 to accomodate the energy // identity and energy solution snapshot. // We add the FOM snapshots in the bases so that they are used // in deriving the EQP rule. - // Also, the change is made here so that the right basis + // Also, the change is made here so that the right basis // dimensions are written in the window parameters file in the // caller's scope. - cutoff[1] += 1; - cutoff[2] += 2; - } - + cutoff[1] += 1; + cutoff[2] += 2; + } + MPI_Bcast(cutoff.GetData(), cutoff.Size(), MPI_INT, 0, MPI_COMM_WORLD); - const CAROM::Matrix *snapX = generator_X->getSnapshotMatrix(); - const CAROM::Matrix *snapV = generator_V->getSnapshotMatrix(); - const CAROM::Matrix *snapE = generator_E->getSnapshotMatrix(); + const CAROM::Matrix *basisV = generator_V->getSpatialBasis(); + const CAROM::Matrix *basisE = generator_E->getSpatialBasis(); - const CAROM::Matrix *basisV = generator_V->getSpatialBasis(); - const CAROM::Matrix *basisE = generator_E->getSpatialBasis(); + // Form the reduced bases. + CAROM::Matrix *tBasisV = new CAROM::Matrix(tH1size, cutoff[1], true); + CAROM::Matrix *tBasisE = new CAROM::Matrix(tL2size, cutoff[2], true); - // Form the reduced bases. - CAROM::Matrix *tBasisV = new CAROM::Matrix(tH1size, cutoff[1], true); - CAROM::Matrix *tBasisE = new CAROM::Matrix(tL2size, cutoff[2], true); - if (input.window == 0) { // Get the first cutoff[1] columns of basisV @@ -1143,7 +1143,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, unitE = 1.0; for (int i = 0; i < tL2size; i++) (*tBasisE)(i, cutoff[2] - 1) = unitE[i]; - + tBasisE->orthogonalize_last(); } else if (input.window > 0) @@ -1159,10 +1159,10 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, // previous window's last snapshot, since we are not using // offset vectors. for (int i = 0; i < tH1size; i++) - (*tBasisV)(i, cutoff[1] - 1) = (*snapV)(i, 0); - + (*tBasisV)(i, cutoff[1] - 1) = Vsnap0(i, 0); + tBasisV->orthogonalize_last(); - + // Get the first cutoff[2] - 2 columns of basisE for (int i = 0; i < tL2size; i++) for (int j = 0; j < cutoff[2] - 2; j++) @@ -1174,25 +1174,25 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, unitE = 1.0; for (int i = 0; i < tL2size; i++) (*tBasisE)(i, cutoff[2] - 2) = unitE[i]; - + tBasisE->orthogonalize_last(cutoff[2] - 1); - + // Add the first E snapshot as the last E basis column // and reorthonormalize. // The current window's first snapshot is the same as the // previous window's last snapshot, since we are not using // offset vectors. for (int i = 0; i < tL2size; i++) - (*tBasisE)(i, cutoff[2] - 1) = (*snapE)(i, 0); - + (*tBasisE)(i, cutoff[2] - 1) = Esnap0(i, 0); + tBasisE->orthogonalize_last(); } - SetupEQP_Force(snapX, snapV, snapE, tBasisV, tBasisE, input, sol); + SetupEQP_Force(&Xsnap0, &Vsnap0, &Esnap0, + tBasisV, tBasisE, input, sol); - // delete snapX, snapV, snapE; - delete tBasisV, tBasisE; - } + delete tBasisV, tBasisE; + } delete generator_X; delete generator_V; From 90b38be2788ab8ef1e133ab5971f88a084fdafa4 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Thu, 7 Dec 2023 08:17:07 -0800 Subject: [PATCH 31/39] basis orthonormality tests --- rom/laghos_rom.cpp | 48 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 47 insertions(+), 1 deletion(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 8dafd4ee..a43cbc6f 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -2299,7 +2299,7 @@ void ROM_Basis::ReadSolutionBases(const int window) { if (window == 0) { - // For the first window, read all rdimv basis vectors normally. + // For the first window, read all rdimv basis vectors normally. basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + std::to_string(window) + basisIdentifier, tH1size, rdimv); } @@ -2365,6 +2365,18 @@ void ROM_Basis::ReadSolutionBases(const int window) (*basisE)(i, tmp_rdime) = unitE[i]; } basisE->orthogonalize_last(); + + std::cout << "Ortho test 1" << "\n"; + Vector test_ortho(rdime); + test_ortho = 0.0; + + for (int j = 0; j < rdime; j++) + { + for (int i = 0; i < tL2size; i++) + test_ortho[j] += (*basisE)(i, j) * (*basisE)(i, rdime - 1); + + std::cout << test_ortho[j] << "\n"; + } } else { @@ -2591,6 +2603,18 @@ void ROM_Basis::AddLastCol_V(Vector const& f) (*basisV)(i, rdimv-1) = mfH1[i]; basisV->orthogonalize_last(); + + std::cout << "Ortho test 2" << "\n"; + Vector test_ortho(rdimv); + test_ortho = 0.0; + + for (int j = 0; j < rdimv; j++) + { + for (int i = 0; i < tH1size; i++) + test_ortho[j] += (*basisV)(i, j) * (*basisV)(i, rdimv - 1); + + std::cout << test_ortho[j] << "\n"; + } } // f is a full vector, not a true vector @@ -2607,6 +2631,28 @@ void ROM_Basis::AddLastCol_E(Vector const& f) (*basisE)(i, rdime-1) = mfL2[i]; basisE->orthogonalize_last(); + + std::cout << "Ortho test 3a" << "\n"; + Vector test_ortho(rdime); + test_ortho = 0.0; + + for (int j = 0; j < rdime - 1; j++) + { + for (int i = 0; i < tL2size; i++) + test_ortho[j] += (*basisE)(i, j) * (*basisE)(i, rdime - 2); + + std::cout << test_ortho[j] << "\n"; + } + + std::cout << "Ortho test 3b" << "\n"; + test_ortho = 0.0; + for (int j = 0; j < rdime; j++) + { + for (int i = 0; i < tL2size; i++) + test_ortho[j] += (*basisE)(i, j) * (*basisE)(i, rdime - 1); + + std::cout << test_ortho[j] << "\n"; + } } // f is a full vector, not a true vector From ea2cdd57ab93cc855e2d94bd2b193b12b734eb6e Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Wed, 20 Dec 2023 11:43:15 -0800 Subject: [PATCH 32/39] Dylan's quadrature data bug fix --- rom/laghos.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 4afabde9..b9ffb206 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -2062,6 +2062,8 @@ int main(int argc, char *argv[]) romOper[romOptions.window-1]->PostprocessHyperreduction(romS); } + if (fom_data) oper->ResetQuadratureData(); + int rdimxprev = romOptions.dimX; int rdimvprev = romOptions.dimV; int rdimeprev = romOptions.dimE; @@ -2095,6 +2097,7 @@ int main(int argc, char *argv[]) onlinePreprocessTimer.Stop(); timeLoopTimer.Start(); } + if (romOptions.use_sample_mesh) { basis[romOptions.window]->ProjectFromPreviousWindow(romOptions, romS, romOptions.window, rdimxprev, rdimvprev, rdimeprev); From 317ca533d98b8752f3e85c57849451f0a8461c43 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Mon, 22 Jan 2024 10:06:57 -0800 Subject: [PATCH 33/39] use double-pass GS --- rom/laghos_rom.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index a43cbc6f..a16bcf72 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -1144,7 +1144,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, for (int i = 0; i < tL2size; i++) (*tBasisE)(i, cutoff[2] - 1) = unitE[i]; - tBasisE->orthogonalize_last(); + tBasisE->orthogonalize_last(-1, true); } else if (input.window > 0) { @@ -1161,7 +1161,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, for (int i = 0; i < tH1size; i++) (*tBasisV)(i, cutoff[1] - 1) = Vsnap0(i, 0); - tBasisV->orthogonalize_last(); + tBasisV->orthogonalize_last(-1, true); // Get the first cutoff[2] - 2 columns of basisE for (int i = 0; i < tL2size; i++) @@ -1175,7 +1175,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, for (int i = 0; i < tL2size; i++) (*tBasisE)(i, cutoff[2] - 2) = unitE[i]; - tBasisE->orthogonalize_last(cutoff[2] - 1); + tBasisE->orthogonalize_last(cutoff[2] - 1, true); // Add the first E snapshot as the last E basis column // and reorthonormalize. @@ -1185,7 +1185,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, for (int i = 0; i < tL2size; i++) (*tBasisE)(i, cutoff[2] - 1) = Esnap0(i, 0); - tBasisE->orthogonalize_last(); + tBasisE->orthogonalize_last(-1, true); } SetupEQP_Force(&Xsnap0, &Vsnap0, &Esnap0, @@ -2364,7 +2364,7 @@ void ROM_Basis::ReadSolutionBases(const int window) { (*basisE)(i, tmp_rdime) = unitE[i]; } - basisE->orthogonalize_last(); + basisE->orthogonalize_last(-1, true); std::cout << "Ortho test 1" << "\n"; Vector test_ortho(rdime); @@ -2404,7 +2404,7 @@ void ROM_Basis::ReadSolutionBases(const int window) { (*basisE)(i, tmp_rdime) = unitE[i]; } - basisE->orthogonalize_last(rdime-1); + basisE->orthogonalize_last(rdime-1, true); // The addition of the lifted energy solution vector as the last // column vector takes place later in the main driver "laghos.cpp" @@ -2602,7 +2602,7 @@ void ROM_Basis::AddLastCol_V(Vector const& f) for (int i = 0; i < tH1size; i++) (*basisV)(i, rdimv-1) = mfH1[i]; - basisV->orthogonalize_last(); + basisV->orthogonalize_last(-1, true); std::cout << "Ortho test 2" << "\n"; Vector test_ortho(rdimv); @@ -2630,7 +2630,7 @@ void ROM_Basis::AddLastCol_E(Vector const& f) for (int i = 0; i < tL2size; i++) (*basisE)(i, rdime-1) = mfL2[i]; - basisE->orthogonalize_last(); + basisE->orthogonalize_last(-1, true); std::cout << "Ortho test 3a" << "\n"; Vector test_ortho(rdime); From 332d59cfec1241b9913007d383379b12da3557cf Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Tue, 6 Feb 2024 14:22:31 -0800 Subject: [PATCH 34/39] perform LQ decomposition of NNLS matrix --- rom/laghos_rom.cpp | 138 +++++++++++++++++++++++++-------------------- 1 file changed, 77 insertions(+), 61 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index a16bcf72..f9fbc0b6 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -421,23 +421,83 @@ void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, { CAROM::NNLSSolver nnls(nnls_tol, 0, maxNNLSnnz, 2); - CAROM::Vector rhs_ub(Gt.numColumns(), false); - // G.mult(w, rhs_ub); // rhs = Gw - // rhs = Gw. Note that by using Gt and multTranspose, we do parallel communication. - Gt.transposeMult(w, rhs_ub); + // rhs = Gw. + // Note that by using Gt and multTranspose, we do parallel communication. + CAROM::Vector rhs_Gw(Gt.numColumns(), false); + Gt.transposeMult(w, rhs_Gw); + + if (true) + { + // Compute Q^T of the LQ factorization of G. + CAROM::Matrix* Qt_ptr; + Qt_ptr = Gt.qr_factorize(); + + CAROM::Matrix Qt(Qt_ptr->getData(), Qt_ptr->numRows(), + Qt_ptr->numColumns(), Qt_ptr->distributed(), + false); + + // Compute L of the factorization; L is lower triangular. + // G = L * Q --> L = G * Q^T + CAROM::Matrix L(Qt.numColumns(), Qt.numColumns(), false); + Gt.transposeMult(Qt, L); + + // Compute the RHS vector. + CAROM::Vector rhs_ub(Qt.numColumns(), false); + Qt.transposeMult(w, rhs_ub); + + // Compute the new RHS tolerance values. + const double delta = 1.0e-11; + CAROM::Vector delta_new(rhs_ub.dim(), false); + for (int i = 0; i < delta_new.dim(); ++i) + { + double denominator = (i + 1) * std::abs(L.item(i, i)); + if (std::abs(denominator) < delta) + delta_new(i) = 1.0; + else + delta_new(i) = delta / denominator; - CAROM::Vector rhs_lb(rhs_ub); - CAROM::Vector rhs_Gw(rhs_ub); + for (int j = i + 1; j < delta_new.dim(); ++j) + { + denominator = (j + 1) * std::abs(L.item(j, i)); - const double delta = 1.0e-11; - for (int i=0; iorthogonalize_last(-1, true); - - std::cout << "Ortho test 1" << "\n"; - Vector test_ortho(rdime); - test_ortho = 0.0; - - for (int j = 0; j < rdime; j++) - { - for (int i = 0; i < tL2size; i++) - test_ortho[j] += (*basisE)(i, j) * (*basisE)(i, rdime - 1); - - std::cout << test_ortho[j] << "\n"; - } } else { @@ -2603,18 +2653,6 @@ void ROM_Basis::AddLastCol_V(Vector const& f) (*basisV)(i, rdimv-1) = mfH1[i]; basisV->orthogonalize_last(-1, true); - - std::cout << "Ortho test 2" << "\n"; - Vector test_ortho(rdimv); - test_ortho = 0.0; - - for (int j = 0; j < rdimv; j++) - { - for (int i = 0; i < tH1size; i++) - test_ortho[j] += (*basisV)(i, j) * (*basisV)(i, rdimv - 1); - - std::cout << test_ortho[j] << "\n"; - } } // f is a full vector, not a true vector @@ -2631,28 +2669,6 @@ void ROM_Basis::AddLastCol_E(Vector const& f) (*basisE)(i, rdime-1) = mfL2[i]; basisE->orthogonalize_last(-1, true); - - std::cout << "Ortho test 3a" << "\n"; - Vector test_ortho(rdime); - test_ortho = 0.0; - - for (int j = 0; j < rdime - 1; j++) - { - for (int i = 0; i < tL2size; i++) - test_ortho[j] += (*basisE)(i, j) * (*basisE)(i, rdime - 2); - - std::cout << test_ortho[j] << "\n"; - } - - std::cout << "Ortho test 3b" << "\n"; - test_ortho = 0.0; - for (int j = 0; j < rdime; j++) - { - for (int i = 0; i < tL2size; i++) - test_ortho[j] += (*basisE)(i, j) * (*basisE)(i, rdime - 1); - - std::cout << test_ortho[j] << "\n"; - } } // f is a full vector, not a true vector From 42f9b0c27d80ee36fd9f4bee5a79052ad597b575 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Sun, 3 Mar 2024 18:13:20 -0800 Subject: [PATCH 35/39] add offset vectors to basis matrices --- rom/laghos.cpp | 47 ++-- rom/laghos_rom.cpp | 561 ++++++++++++++++++++++++++----------------- rom/laghos_rom.hpp | 19 +- rom/laghos_utils.cpp | 7 +- 4 files changed, 372 insertions(+), 262 deletions(-) diff --git a/rom/laghos.cpp b/rom/laghos.cpp index b9ffb206..2abeaa43 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -1262,13 +1262,14 @@ int main(int argc, char *argv[]) else { // Construct the ROM basis and operator only for the first - // window; the ones for the next windows will be constructed - // at the time of each window change, because the current - // vector S is needed for their construction. - // In particular, forming the bases and reduced inverse mass - // matrices requires knoweldege of vector S. + // window. + // Those for the next windows will be constructed at the + // time of each window change, because the current solution + // vector S must be known. SetWindowParameters(twparam, romOptions); - basis[0] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV, ×teps); + basis[0] = new ROM_Basis(romOptions, MPI_COMM_WORLD, + sFactorX, sFactorV, ×teps, S); + if (!romOptions.hyperreduce_prep) { romOper[0] = new ROM_Operator(romOptions, basis[0], rho_coeff, mat_coeff, order_e, source, visc, vort, cfl, p_assembly, @@ -1278,7 +1279,8 @@ int main(int argc, char *argv[]) } else { - basis[0] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV, ×teps); + basis[0] = new ROM_Basis(romOptions, MPI_COMM_WORLD, + sFactorX, sFactorV, ×teps, S); if (!romOptions.hyperreduce_prep) { romOper[0] = new ROM_Operator(romOptions, basis[0], rho_coeff, mat_coeff, order_e, source, visc, vort, cfl, p_assembly, @@ -1475,7 +1477,8 @@ int main(int argc, char *argv[]) SetWindowParameters(twparam, romOptions); } - basis[0] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV); + basis[0] = new ROM_Basis(romOptions, MPI_COMM_WORLD, + sFactorX, sFactorV, nullptr, S); basis[0]->Init(romOptions, *S); if (romOptions.mergeXV) @@ -1537,15 +1540,8 @@ int main(int argc, char *argv[]) basis[romOptions.window-1]->LiftROMtoFOM(romS, *S); delete basis[romOptions.window-1]; - basis[romOptions.window] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV); - - if (romOptions.hyperreductionSamplingType == eqp_energy) - { - // Add the last lifted solution vector as the last - // column in the bases. - basis[romOptions.window]->AddLastCol_V(*S); - basis[romOptions.window]->AddLastCol_E(*S); - } + basis[romOptions.window] = new ROM_Basis(romOptions, + MPI_COMM_WORLD, sFactorX, sFactorV, nullptr, S); basis[romOptions.window]->Init(romOptions, *S); if (romOptions.mergeXV) @@ -2075,19 +2071,10 @@ int main(int argc, char *argv[]) onlinePreprocessTimer.Start(); // Form the ROM basis and operator for the new window. - basis[romOptions.window] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV, ×teps); - - onlinePreprocessTimer.Stop(); - timeLoopTimer.Start(); - - // Add the last lifted solution vector as the last - // column in the bases. - basis[romOptions.window]->AddLastCol_V(*S); - basis[romOptions.window]->AddLastCol_E(*S); - - timeLoopTimer.Stop(); - onlinePreprocessTimer.Start(); - + basis[romOptions.window] = new ROM_Basis(romOptions, + MPI_COMM_WORLD, sFactorX, sFactorV, + ×teps, S); + if (!romOptions.hyperreduce_prep) { romOper[romOptions.window] = new ROM_Operator(romOptions, basis[romOptions.window], rho_coeff, mat_coeff, diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index f9fbc0b6..b18e19fa 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -441,6 +441,22 @@ void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, CAROM::Matrix L(Qt.numColumns(), Qt.numColumns(), false); Gt.transposeMult(Qt, L); + // Check for nearly linearly dependent Q rows. + // This is achieved by checking the magnitude of the diagonal + // values of L. + std::vector row_ind; + for (int i = 0; i < L.numRows(); ++i) + { + if (std::abs(L.item(i, i)) < 1e-12) + { + row_ind.push_back(i); + std::cout << i << ", "; + } + } + std::cout << "\n"; + std::cout << "Found " << row_ind.size() << " / " << + L.numRows() << " nearly linearly dependent constraints.\n"; + // Compute the RHS vector. CAROM::Vector rhs_ub(Qt.numColumns(), false); Qt.transposeMult(w, rhs_ub); @@ -1017,7 +1033,7 @@ void ROM_Sampler::SetupEQP_Force(const CAROM::Matrix* snapX, const CAROM::Matrix* basisV, const CAROM::Matrix* basisE, ROM_Options const& input, - Vector const& sol) + BlockVector const& sol) { MFEM_VERIFY(basisV->numRows() == input.H1FESpace->GetTrueVSize(), ""); MFEM_VERIFY(basisE->numRows() == input.L2FESpace->GetTrueVSize(), ""); @@ -1043,7 +1059,7 @@ void ROM_Sampler::SetupEQP_Force(const CAROM::Matrix* snapX, } void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, - Vector const& sol) + BlockVector const& sol) { CAROM::Matrix Xsnap0(*(generator_X->getSnapshotMatrix())); CAROM::Matrix Vsnap0(*(generator_V->getSnapshotMatrix())); @@ -1156,100 +1172,135 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, if (input.hyperreductionSamplingType == eqp_energy) { - if (rank == 0 && input.window == 0) - { - // If first window, increase the energy basis dimension by 1 - // to accomodate the energy identity. - // No change in the velocity bases dimension. - cutoff[2] += 1; - } - else if (rank == 0 && input.window > 0) + if (rank == 0) { - // If not first window, increase the dimension of the velocity - // basis by 1 to accomodate the velocity solution snapshot, - // and the energy basis dimension by 2 to accomodate the energy - // identity and energy solution snapshot. - // We add the FOM snapshots in the bases so that they are used - // in deriving the EQP rule. - // Also, the change is made here so that the right basis - // dimensions are written in the window parameters file in the - // caller's scope. - cutoff[1] += 1; - cutoff[2] += 2; + // Increase the X basis dimension by 1 to accomodate the + // addition of the lifted X solution in the online stage + // (the X basis is not needed for the EQP rule calculation). + cutoff[0] += 1; + + // Increase the V basis dimension by 1 to accomodate the + // V offset vector. + // The V offset vector is added to the basis only if it + // has nonzero norm. + if ( !(input.window == 0 && initV->norm() < 1e-15) ) + { + cutoff[1] += 1; + } + + // Increase the E basis dimension by 2 to accomodate the + // E identity and E offset vector. + // The E offset vector is added to the basis only if it + // has nonzero norm. + if (input.window == 0 && initE->norm() < 1e-15) + { + cutoff[2] += 1; + } + else + { + cutoff[2] += 2; + } + + // Offsets are used only in the offline stage, to aid in the + // efficient derivation of reduced bases and EQP rule + // construction. + // Once the reduced bases are derived, the offsets are added + // as the last basis column. + // In the online stage, the appropriate lifted solution vectors + // are instead added to the reduced bases and the simulation + // is run without the use of offset vectors. + + // The changes are made here so that the right basis dimensions + // are written in the window parameters file in the caller's + // scope. } MPI_Bcast(cutoff.GetData(), cutoff.Size(), MPI_INT, 0, MPI_COMM_WORLD); - const CAROM::Matrix *basisV = generator_V->getSpatialBasis(); - const CAROM::Matrix *basisE = generator_E->getSpatialBasis(); + // Get the bases resulting from the SVDs. + CAROM::Matrix const* basisV = generator_V->getSpatialBasis(); + CAROM::Matrix const* basisE = generator_E->getSpatialBasis(); - // Form the reduced bases. - CAROM::Matrix *tBasisV = new CAROM::Matrix(tH1size, cutoff[1], true); - CAROM::Matrix *tBasisE = new CAROM::Matrix(tL2size, cutoff[2], true); + CAROM::Matrix* tBasisV = nullptr; + CAROM::Matrix* tBasisE = nullptr; - if (input.window == 0) + // V basis construction. + if (input.window == 0 && initV->norm() < 1e-15) { - // Get the first cutoff[1] columns of basisV - for (int i = 0; i < tH1size; i++) - for (int j = 0; j < cutoff[1]; j++) + // Get the first cutoff[1] columns of basisV. + tBasisV = basisV->getFirstNColumns(cutoff[1]); + } + else + { + tBasisV = new CAROM::Matrix(tH1size, cutoff[1], true); + + // Get the first cutoff[1] - 1 columns of basisV. + for (int i = 0; i < tH1size; ++i) + { + for (int j = 0; j < cutoff[1] - 1; ++j) + { (*tBasisV)(i, j) = (*basisV)(i, j); + } + } + // Add the V offset vector as the last V basis column + // and reorthonormalize. + for (int i = 0; i < tH1size; ++i) + { + (*tBasisV)(i, cutoff[1] - 1) = (*initV)(i); + } + tBasisV->orthogonalize_last(-1, true); + } - // Get the first cutoff[2] - 1 columns of basisE - for (int i = 0; i < tL2size; i++) - for (int j = 0; j < cutoff[2] - 1; j++) - (*tBasisE)(i, j) = (*basisE)(i, j); + // E basis construction. + if (input.window == 0 && initE->norm() < 1e-15) + { + tBasisE = new CAROM::Matrix(tL2size, cutoff[2], true); + // Get the first cutoff[2] - 1 columns of basisE. + for (int i = 0; i < tL2size; ++i) + { + for (int j = 0; j < cutoff[2] - 1; ++j) + { + (*tBasisE)(i, j) = (*basisE)(i, j); + } + } // Add the energy identity as the last E basis column // and reorthonormalize. - Vector unitE(tL2size); - unitE = 1.0; - for (int i = 0; i < tL2size; i++) - (*tBasisE)(i, cutoff[2] - 1) = unitE[i]; - + for (int i = 0; i < tL2size; ++i) + { + (*tBasisE)(i, cutoff[2] - 1) = 1.0; + } tBasisE->orthogonalize_last(-1, true); } - else if (input.window > 0) + else { - // Get the first cutoff[1] - 1 columns of basisV - for (int i = 0; i < tH1size; i++) - for (int j = 0; j < cutoff[1] - 1; j++) - (*tBasisV)(i, j) = (*basisV)(i, j); - - // Add the first V snapshot as the last V basis column - // and reorthonormalize. - // The current window's first snapshot is the same as the - // previous window's last snapshot, since we are not using - // offset vectors. - for (int i = 0; i < tH1size; i++) - (*tBasisV)(i, cutoff[1] - 1) = Vsnap0(i, 0); - - tBasisV->orthogonalize_last(-1, true); + tBasisE = new CAROM::Matrix(tL2size, cutoff[2], true); - // Get the first cutoff[2] - 2 columns of basisE - for (int i = 0; i < tL2size; i++) - for (int j = 0; j < cutoff[2] - 2; j++) + // Get the first cutoff[2] - 2 columns of basisE. + for (int i = 0; i < tL2size; ++i) + { + for (int j = 0; j < cutoff[2] - 2; ++j) + { (*tBasisE)(i, j) = (*basisE)(i, j); - + } + } // Add the energy identity as the penultimate E basis column // and reorthonormalize. - Vector unitE(tL2size); - unitE = 1.0; - for (int i = 0; i < tL2size; i++) - (*tBasisE)(i, cutoff[2] - 2) = unitE[i]; - + for (int i = 0; i < tL2size; ++i) + { + (*tBasisE)(i, cutoff[2] - 2) = 1.0; + } tBasisE->orthogonalize_last(cutoff[2] - 1, true); - // Add the first E snapshot as the last E basis column + // Add the E offset vector as the last E basis column // and reorthonormalize. - // The current window's first snapshot is the same as the - // previous window's last snapshot, since we are not using - // offset vectors. - for (int i = 0; i < tL2size; i++) - (*tBasisE)(i, cutoff[2] - 1) = Esnap0(i, 0); - + for (int i = 0; i < tL2size; ++i) + { + (*tBasisE)(i, cutoff[2] - 1) = (*initE)(i); + } tBasisE->orthogonalize_last(-1, true); } - + SetupEQP_Force(&Xsnap0, &Vsnap0, &Esnap0, tBasisV, tBasisE, input, sol); @@ -1342,8 +1393,10 @@ CAROM::Matrix* MultBasisROM(const int rank, const std::string filename, const in return S; } -ROM_Basis::ROM_Basis(ROM_Options const& input, MPI_Comm comm_, const double sFactorX, const double sFactorV, - const std::vector *timesteps) +ROM_Basis::ROM_Basis(ROM_Options const& input, MPI_Comm comm_, + const double sFactorX, const double sFactorV, + const std::vector *timesteps, + BlockVector const* S) : comm(comm_), rdimx(input.dimX), rdimv(input.dimV), rdime(input.dimE), rdimfv(input.dimFv), rdimfe(input.dimFe), numSamplesX(input.sampX), numSamplesV(input.sampV), numSamplesE(input.sampE), numTimeSamplesV(input.tsampV), numTimeSamplesE(input.tsampE), @@ -1402,42 +1455,7 @@ ROM_Basis::ROM_Basis(ROM_Options const& input, MPI_Comm comm_, const double sFac mfH1.SetSize(tH1size); mfL2.SetSize(tL2size); - // This code block is used to set the sizes of the basis matrices - // correctly, before the reading of the actual data takes place - // within ReadSolutionBases following right after. - // In this way, once the basis data has been read, the matrices - // can be extended by adding the required column vectors. - if (hyperreductionSamplingType == eqp_energy) - { - if (input.window == 0) - { - // For the first window, only the energy basis need be extended. - - // The energy basis is extended by adding the energy identity. - // Dimension rdime has already been increased by 1 in the - // offline stage - basisE = new CAROM::Matrix(tL2size, rdime, true); - } - else - { - // For any window other than the first, both energy and - // velocity bases need be extended. - - // The velocity basis is extended by adding the current lifted - // velocity solution vector. - // Dimension rdimv has already been increased by 1 in the - // offline stage. - basisV = new CAROM::Matrix(tH1size, rdimv, true); - - // The energy basis is extended by adding the energy identity - // and the current lifted energy solution vector. - // Dimension rdime has already been increased by 2 in the - // offline stage. - basisE = new CAROM::Matrix(tL2size, rdime, true); - } - } - - ReadSolutionBases(input.window); + ReadSolutionBases(input.window, S); if (spaceTime) { @@ -2354,41 +2372,62 @@ int ROM_Basis::SolutionSizeFOM() const return (2*H1size) + L2size; // full size, not true DOF size } -void ROM_Basis::ReadSolutionBases(const int window) +void ROM_Basis::ReadSolutionBases(const int window, BlockVector const* S) { // Velocity basis formation. if (hyperreductionSamplingType == eqp_energy) { - if (window == 0) - { - // For the first window, read all rdimv basis vectors normally. - basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + - std::to_string(window) + basisIdentifier, tH1size, rdimv); - } - else - { - // For any window other than the first, read the first rdimv-1 - // basis vectors, since rdimv has been increased by 1 to accomodate - // the addition of the lifted velocity solution vector. - int tmp_rdimv = rdimv - 1; - CAROM::Matrix *tmp_basisV = 0; - - tmp_basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + - std::to_string(window) + basisIdentifier, tH1size, tmp_rdimv); + if (window == 0) + { + // TODO: use S to check the initial condition norm directly, + // rather than reading the stored file. + initV = new CAROM::Vector(tH1size, true); + std::string path_init = testing_parameter_basename + + "/ROMoffset" + basisIdentifier + "/init"; + cout << "Reading: " << path_init << endl; + initV->read(path_init + "V0"); - for (int i = 0; i < tH1size; i++) - { - for (int j = 0; j < tmp_rdimv; j++) - (*basisV)(i, j) = (*tmp_basisV)(i, j); - } - delete tmp_basisV; - - // The addition of the lifted velocity solution vector as the last - // column vector takes place later in the main driver "laghos.cpp" - // at the time of the window change, because the current solution - // vector S must be known. - } - } + if (initV->norm() < 1e-15) + { + extendFirstWindowV = false; // true by default + } + //delete initV; // gets deleted by ROM_Basis destructor + } + if (window == 0 && !extendFirstWindowV) + { + // The V basis is not extended. + basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + + std::to_string(window) + basisIdentifier, + tH1size, rdimv); + } + else + { + // The V basis is extended by adding the current lifted + // V solution vector. + // Dimension rdimv has already been set accordingly + // in the offline stage. + basisV = new CAROM::Matrix(tH1size, rdimv, true); + + // Read the first (rdimv - 1) basis vectors. + int tmp_rdimv = rdimv - 1; + CAROM::Matrix* tmp_basisV = nullptr; + + tmp_basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + + std::to_string(window) + basisIdentifier, + tH1size, tmp_rdimv); + + // Add them to the V basis. + for (int i = 0; i < tH1size; i++) + { + for (int j = 0; j < tmp_rdimv; j++) + (*basisV)(i, j) = (*tmp_basisV)(i, j); + } + delete tmp_basisV; + + // Add the lifted V solution as the last basis vector. + AddLastCol_V(*S); + } + } else { if (!useVX) @@ -2399,80 +2438,135 @@ void ROM_Basis::ReadSolutionBases(const int window) } // Energy basis formation. - if (hyperreductionSamplingType == eqp_energy) - { - if (window == 0) - { - // For the first window, read the first rdime-1 basis vectors, - // since rdime has been increased by 1 to accomodate the addition - // of the energy identity. - int tmp_rdime = rdime - 1; - CAROM::Matrix *tmp_basisE = 0; - - tmp_basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + - std::to_string(window) + basisIdentifier, tL2size, tmp_rdime); - - for (int i = 0; i < tL2size; i++) - { - for (int j = 0; j < tmp_rdime; j++) - (*basisE)(i, j) = (*tmp_basisE)(i, j); - } - delete tmp_basisE; - - // Add the energy identity as the last column vector. - Vector unitE(tL2size); - unitE = 1.0; - for (int i = 0; i < tL2size; i++) - { - (*basisE)(i, tmp_rdime) = unitE[i]; - } - basisE->orthogonalize_last(-1, true); - } - else - { - // For any window other than the first, read the first rdime-2 - // basis vectors, since rdime has been increased by 2 to accomodate - // the addition of the energy identity and the lifted energy - // solution vector. - int tmp_rdime = rdime - 2; - CAROM::Matrix *tmp_basisE = 0; - - tmp_basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + - std::to_string(window) + basisIdentifier, tL2size, tmp_rdime); - - for (int i = 0; i < tL2size; i++) - { - for (int j = 0; j < tmp_rdime; j++) - (*basisE)(i, j) = (*tmp_basisE)(i, j); - } - delete tmp_basisE; - - // Add the energy identity as the penultimate column vector. - Vector unitE(tL2size); - unitE = 1.0; - for (int i = 0; i < tL2size; i++) - { - (*basisE)(i, tmp_rdime) = unitE[i]; - } - basisE->orthogonalize_last(rdime-1, true); - - // The addition of the lifted energy solution vector as the last - // column vector takes place later in the main driver "laghos.cpp" - // at the time of the window change, because the current solution - // vector S must be known. - } - } + if (hyperreductionSamplingType == eqp_energy) + { + if (window == 0) + { + // TODO: use S to check the initial condition norm directly, + // rather than reading the stored file. + initE = new CAROM::Vector(tL2size, true); + std::string path_init = testing_parameter_basename + + "/ROMoffset" + basisIdentifier + "/init"; + cout << "Reading: " << path_init << endl; + initE->read(path_init + "E0"); + + if (initE->norm() < 1e-15) + { + extendFirstWindowE = false; // true by default + } + //delete initE; // gets deleted by ROM_Basis destructor + } + if (window == 0 && !extendFirstWindowE) + { + // The E basis is extended by adding the E identity. + // Dimension rdime has already been set accordingly + // in the offline stage. + basisE = new CAROM::Matrix(tL2size, rdime, true); + + // Read the first (rdime - 1) basis vectors. + int tmp_rdime = rdime - 1; + CAROM::Matrix* tmp_basisE = nullptr; + + tmp_basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + + std::to_string(window) + basisIdentifier, tL2size, + tmp_rdime); + + // Add them to the E basis. + for (int i = 0; i < tL2size; i++) + { + for (int j = 0; j < tmp_rdime; j++) + (*basisE)(i, j) = (*tmp_basisE)(i, j); + } + delete tmp_basisE; + + // Add the E identity as the last column vector. + for (int i = 0; i < tL2size; i++) + { + (*basisE)(i, tmp_rdime) = 1.0; + } + basisE->orthogonalize_last(-1, true); + } + else + { + // The E basis is extended by adding the E identity and the + // current lifted E solution vector. + // Dimension rdime has already been set accordingly + // in the offline stage. + basisE = new CAROM::Matrix(tL2size, rdime, true); + + // Read the first (rdime - 2) basis vectors. + int tmp_rdime = rdime - 2; + CAROM::Matrix* tmp_basisE = nullptr; + + tmp_basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + + std::to_string(window) + basisIdentifier, tL2size, + tmp_rdime); + + // Add them to the E basis. + for (int i = 0; i < tL2size; i++) + { + for (int j = 0; j < tmp_rdime; j++) + (*basisE)(i, j) = (*tmp_basisE)(i, j); + } + delete tmp_basisE; + + // Add the E identity as the penultimate column vector. + for (int i = 0; i < tL2size; i++) + { + (*basisE)(i, tmp_rdime) = 1.0; + } + basisE->orthogonalize_last(rdime-1, true); + + // Add the lifted E solution as the last basis vector. + AddLastCol_E(*S); + } + } else { basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + std::to_string(window) + basisIdentifier, tL2size, rdime); } + // Position basis formation. + if (hyperreductionSamplingType == eqp_energy) + { + // The X basis is extended by adding the current lifted + // X solution vector. + // Dimension rdimx has already been set accordingly + // in the offline stage. + basisX = new CAROM::Matrix(tH1size, rdimx, true); + + // Read the first (rdimx - 1) basis vectors. + int tmp_rdimx = rdimx - 1; + CAROM::Matrix* tmp_basisX = nullptr; + + tmp_basisX = ReadBasisROM(rank, basename + "/" + ROMBasisName::X + + std::to_string(window) + basisIdentifier, tH1size, tmp_rdimx); - if (useXV) - basisX = basisV; + for (int i = 0; i < tH1size; i++) + { + for (int j = 0; j < tmp_rdimx; j++) + { + (*basisX)(i, j) = (*tmp_basisX)(i, j); + } + } + delete tmp_basisX; + + // Add the lifted X solution as the last basis vector. + AddLastCol_X(*S); + } else - basisX = ReadBasisROM(rank, basename + "/" + ROMBasisName::X + std::to_string(window) + basisIdentifier, tH1size, rdimx); + { + if (useXV) + { + basisX = basisV; + } + else + { + basisX = ReadBasisROM(rank, basename + "/" + ROMBasisName::X + + std::to_string(window) + basisIdentifier, tH1size, rdimx); + } + } if (useVX) basisV = basisX; @@ -2640,7 +2734,23 @@ void ROM_Basis::ProjectFOMtoROM_V(Vector const& f, Vector & r, const bool timeDe } // f is a full vector, not a true vector -void ROM_Basis::AddLastCol_V(Vector const& f) +void ROM_Basis::AddLastCol_X(BlockVector const& f) +{ + MFEM_VERIFY(f.Size() == (2*H1size) + L2size, ""); + + for (int i = 0; i < H1size; ++i) + (*gfH1)(i) = f[i]; + + gfH1->GetTrueDofs(mfH1); + + for (int i = 0; i < tH1size; ++i) + (*basisX)(i, rdimx-1) = mfH1[i]; + + basisX->orthogonalize_last(-1, true); +} + +// f is a full vector, not a true vector +void ROM_Basis::AddLastCol_V(BlockVector const& f) { MFEM_VERIFY(f.Size() == (2*H1size) + L2size, ""); @@ -2656,7 +2766,7 @@ void ROM_Basis::AddLastCol_V(Vector const& f) } // f is a full vector, not a true vector -void ROM_Basis::AddLastCol_E(Vector const& f) +void ROM_Basis::AddLastCol_E(BlockVector const& f) { MFEM_VERIFY(f.Size() == (2*H1size) + L2size, ""); @@ -4984,13 +5094,16 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res, } } // Loop (i) over EQP points - if (energy_conserve) - { - // Multiply by the reduced Mv inverse - Vector res_tmp = res; - - invMvROM.Mult(res_tmp, res); - } + if (energy_conserve) + { + // Multiply by the reduced Mv inverse + Vector res_tmp(res.Size()); + for (int i = 0; i < res.Size(); ++i) + { + res_tmp[i] = res[i]; + } + invMvROM.Mult(res_tmp, res); + } MPI_Allreduce(MPI_IN_PLACE, res.GetData(), res.Size(), MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); @@ -5171,16 +5284,16 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res, } } // Loop (i) over EQP points - if (energy_conserve) - { - // Multiply by the reduced Me inverse - - // TODO: does that copy the data of res to res_tmp, or is it just - // a reference? - Vector res_tmp = res; - - invMeROM.Mult(res_tmp, res); - } + if (energy_conserve) + { + // Multiply by the reduced Me inverse + Vector res_tmp(res.Size()); + for (int i = 0; i < res.Size(); ++i) + { + res_tmp[i] = res[i]; + } + invMeROM.Mult(res_tmp, res); + } MPI_Allreduce(MPI_IN_PLACE, res.GetData(), res.Size(), MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); diff --git a/rom/laghos_rom.hpp b/rom/laghos_rom.hpp index 9c9b1c7c..55f95123 100644 --- a/rom/laghos_rom.hpp +++ b/rom/laghos_rom.hpp @@ -632,7 +632,7 @@ class ROM_Sampler void SampleSolution(const double t, const double dt, const double pd, Vector const& S); - void Finalize(Array &cutoff, ROM_Options& input, Vector const& sol); + void Finalize(Array &cutoff, ROM_Options& input, BlockVector const& sol); int MaxNumSamples() { @@ -827,7 +827,7 @@ class ROM_Sampler void SetupEQP_Force(const CAROM::Matrix* snapX, const CAROM::Matrix* snapV, const CAROM::Matrix* snapE, const CAROM::Matrix* basisV, const CAROM::Matrix* basisE, ROM_Options const& input, - Vector const& sol); + BlockVector const& sol); void SetupEQP_Force_Eq(const CAROM::Matrix* snapX, const CAROM::Matrix* snapV, const CAROM::Matrix* snapE, const CAROM::Matrix* basisV, const CAROM::Matrix* basisE, ROM_Options const& input, @@ -844,13 +844,14 @@ class ROM_Basis public: ROM_Basis(ROM_Options const& input, MPI_Comm comm_, const double sFactorX=1.0, const double sFactorV=1.0, - const std::vector *timesteps=NULL); + const std::vector *timesteps=NULL, + BlockVector const* S = nullptr); ~ROM_Basis(); void Init(ROM_Options const& input, Vector const& S); - void ReadSolutionBases(const int window); + void ReadSolutionBases(const int window, BlockVector const* S = nullptr); void ReadTemporalBases(const int window); void ProjectFOMtoROM(Vector const& f, Vector & r, @@ -859,8 +860,9 @@ class ROM_Basis void ProjectFOMtoROM_V(Vector const& f, Vector & r, const bool timeDerivative=false); - void AddLastCol_V(Vector const& f); - void AddLastCol_E(Vector const& f); + void AddLastCol_X(BlockVector const& f); + void AddLastCol_V(BlockVector const& f); + void AddLastCol_E(BlockVector const& f); void LiftROMtoFOM(Vector const& r, Vector & f); void LiftROMtoFOM_dVdt(Vector const& r, Vector & f); @@ -992,6 +994,11 @@ class ROM_Basis const bool useVX; // If true, use X-X0 for V. const bool mergeXV; // If true, merge bases for X-X0 and V. + // Those two flags are checked when extending the reduced V & E bases + // in the CEQP hyperreduction case. + bool extendFirstWindowV = true; + bool extendFirstWindowE = true; + int H1size; int L2size; int tH1size; diff --git a/rom/laghos_utils.cpp b/rom/laghos_utils.cpp index 8ac0e232..f7e6c647 100644 --- a/rom/laghos_utils.cpp +++ b/rom/laghos_utils.cpp @@ -178,8 +178,11 @@ void VerifyOfflineParam(int& dim, double& dt, ROM_Options& romOptions, std::getline(opin, line); split_line(line, words); - MFEM_VERIFY(std::stoi(words[2]) == romOptions.useOffset, "-romos option does not match record."); - MFEM_VERIFY(std::stoi(words[3]) == romOptions.offsetType, "-romostype option does not match record."); + if (romOptions.hyperreductionSamplingType != eqp_energy) + { + MFEM_VERIFY(std::stoi(words[2]) == romOptions.useOffset, "-romos option does not match record."); + MFEM_VERIFY(std::stoi(words[3]) == romOptions.offsetType, "-romostype option does not match record."); + } MFEM_VERIFY(std::stoi(words[4]) == romOptions.SNS, "-romsns option does not match record."); if (rom_offline) From e3a7be90624eaba972fb271db718fa8fc714d967 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Mon, 25 Mar 2024 13:59:53 -0700 Subject: [PATCH 36/39] update libROM sampling func. calls --- rom/laghos_rom.cpp | 24 ++++++++++++------------ rom/merge.cpp | 6 +++--- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index b18e19fa..0b521fad 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -82,12 +82,12 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p Xdiff[i] = X[i] - (*initX)(i); } - addSample = generator_X->takeSample(Xdiff.GetData(), t, dt); + addSample = generator_X->takeSample(Xdiff.GetData()); generator_X->computeNextSampleTime(Xdiff.GetData(), dXdt.GetData(), t); } else { - addSample = generator_X->takeSample(X.GetData(), t, dt); + addSample = generator_X->takeSample(X.GetData()); generator_X->computeNextSampleTime(X.GetData(), dXdt.GetData(), t); } @@ -121,12 +121,12 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p Xdiff[i] = V[i] - (*initV)(i); } - addSample = generator_V->takeSample(Xdiff.GetData(), t, dt); + addSample = generator_V->takeSample(Xdiff.GetData()); generator_V->computeNextSampleTime(Xdiff.GetData(), dVdt.GetData(), t); } else { - addSample = generator_V->takeSample(V.GetData(), t, dt); + addSample = generator_V->takeSample(V.GetData()); generator_V->computeNextSampleTime(V.GetData(), dVdt.GetData(), t); } @@ -143,7 +143,7 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p gfH1[i] = dSdt[H1size + i]; // Fv gfH1.GetTrueDofs(Xdiff); - bool addSampleF = generator_Fv->takeSample(Xdiff.GetData(), t, dt); + bool addSampleF = generator_Fv->takeSample(Xdiff.GetData()); if (writeSnapshots && addSampleF) { @@ -170,12 +170,12 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p Ediff[i] = E[i] - (*initE)(i); } - addSample = generator_E->takeSample(Ediff.GetData(), t, dt); + addSample = generator_E->takeSample(Ediff.GetData()); generator_E->computeNextSampleTime(Ediff.GetData(), dEdt.GetData(), t); } else { - addSample = generator_E->takeSample(E.GetData(), t, dt); + addSample = generator_E->takeSample(E.GetData()); generator_E->computeNextSampleTime(E.GetData(), dEdt.GetData(), t); } @@ -191,7 +191,7 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p gfL2[i] = dSdt[(2*H1size) + i]; // Fe gfL2.GetTrueDofs(Ediff); - addSampleF = generator_Fe->takeSample(Ediff.GetData(), t, dt); + addSampleF = generator_Fe->takeSample(Ediff.GetData()); if (writeSnapshots && addSampleF) { @@ -1326,11 +1326,11 @@ CAROM::Matrix* ReadBasisROM(const int rank, const std::string filename, const in CAROM::Matrix *basis; if (dim == -1) { - basis = (CAROM::Matrix*) reader.getSpatialBasis(0.0); + basis = (CAROM::Matrix*) reader.getSpatialBasis(); } else { - basis = (CAROM::Matrix*) reader.getSpatialBasis(0.0, dim); + basis = (CAROM::Matrix*) reader.getSpatialBasis(dim); } MFEM_VERIFY(basis->numRows() == vectorSize, ""); @@ -2588,7 +2588,7 @@ void ROM_Basis::ReadSolutionBases(const int window, BlockVector const* S) ej(j) = 1.0; basisX->mult(ej, CBej); - const bool addSample = generator_XV.takeSample(Bej.GetData(), 0.0, 1.0); // artificial time and timestep + const bool addSample = generator_XV.takeSample(Bej.GetData()); MFEM_VERIFY(addSample, "Sample not added"); ej(j) = 0.0; @@ -2602,7 +2602,7 @@ void ROM_Basis::ReadSolutionBases(const int window, BlockVector const* S) ej(j) = 1.0; basisV->mult(ej, CBej); - const bool addSample = generator_XV.takeSample(Bej.GetData(), 0.0, 1.0); // artificial time and timestep + const bool addSample = generator_XV.takeSample(Bej.GetData()); MFEM_VERIFY(addSample, "Sample not added"); ej(j) = 0.0; diff --git a/rom/merge.cpp b/rom/merge.cpp index 656db5ba..73501e4b 100644 --- a/rom/merge.cpp +++ b/rom/merge.cpp @@ -85,7 +85,7 @@ void MergeSamplingWindow(const int rank, const int first_sv, const double energy int col_ub = std::min(offsetAllWindows[basisWindow+1][paramID+nsets*v]+windowOverlapSamples+1, num_snap); int num_cols = col_ub - col_lb + 1; std::cout << num_cols << " columns read. Columns " << col_lb - 1 << " to " << col_ub - 1 << std::endl; - const CAROM::Matrix* mat = basis_reader->getSnapshotMatrix(0.0, col_lb, col_ub); + const CAROM::Matrix* mat = basis_reader->getSnapshotMatrix(col_lb, col_ub); MFEM_VERIFY(dim == mat->numRows(), "Inconsistent snapshot size"); MFEM_VERIFY(num_cols == mat->numColumns(), "Inconsistent number of snapshots"); @@ -111,7 +111,7 @@ void MergeSamplingWindow(const int rank, const int first_sv, const double energy { tmp[i] = (offsetInit) ? mat->item(i,j) - mat->item(i,0) : mat->item(i,j); } - window_basis_generator->takeSample(tmp.GetData(), 0.0, 1.0); + window_basis_generator->takeSample(tmp.GetData()); } delete mat; @@ -171,7 +171,7 @@ void GetSnapshotDim(const int id, const std::string& basename, const std::string std::string filename = basename + "/param" + std::to_string(id) + "_var" + varName + std::to_string(window) + basisIdentifier + "_snapshot"; CAROM::BasisReader reader(filename); - const CAROM::Matrix *S = reader.getSnapshotMatrix(0.0); + const CAROM::Matrix *S = reader.getSnapshotMatrix(); varDim = S->numRows(); numSnapshots = S->numColumns(); } From 9c8294e7e21ef6a2407ddec098200e0c1cf604b9 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Tue, 21 May 2024 13:37:28 -0700 Subject: [PATCH 37/39] use mass induced Gram-Schmidt --- rom/laghos_rom.cpp | 281 +++++++++++++++++++++++++++---------------- rom/laghos_rom.hpp | 3 +- rom/laghos_utils.cpp | 190 +++++++++++++++++++++++++++++ rom/laghos_utils.hpp | 12 ++ 4 files changed, 384 insertions(+), 102 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index 0b521fad..c9b4e968 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -426,7 +426,10 @@ void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, CAROM::Vector rhs_Gw(Gt.numColumns(), false); Gt.transposeMult(w, rhs_Gw); - if (true) + // TODO: Make it a user input option. + bool const nnls_use_lq = true; + + if (nnls_use_lq) { // Compute Q^T of the LQ factorization of G. CAROM::Matrix* Qt_ptr; @@ -1172,6 +1175,18 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, if (input.hyperreductionSamplingType == eqp_energy) { + // Get the bases resulting from the SVDs. + CAROM::Matrix const* basisV = generator_V->getSpatialBasis(); + CAROM::Matrix const* basisE = generator_E->getSpatialBasis(); + + // TODO: Remove the const qualifier in libROM. + // Then orthonormalize the V and E bases w.r.t. the mass matrices + // and call endSamples() again to rewrite the bases, overwriting + // the previous versions. + // In this way, I will then have to orthonormalize only the new + // vectors I add (including when I read the bases in the + // online stage). + if (rank == 0) { // Increase the X basis dimension by 1 to accomodate the @@ -1214,13 +1229,8 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, // are written in the window parameters file in the caller's // scope. } - MPI_Bcast(cutoff.GetData(), cutoff.Size(), MPI_INT, 0, MPI_COMM_WORLD); - // Get the bases resulting from the SVDs. - CAROM::Matrix const* basisV = generator_V->getSpatialBasis(); - CAROM::Matrix const* basisE = generator_E->getSpatialBasis(); - CAROM::Matrix* tBasisV = nullptr; CAROM::Matrix* tBasisE = nullptr; @@ -1248,7 +1258,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, { (*tBasisV)(i, cutoff[1] - 1) = (*initV)(i); } - tBasisV->orthogonalize_last(-1, true); + //tBasisV->orthogonalize_last(-1, true); } // E basis construction. @@ -1270,7 +1280,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, { (*tBasisE)(i, cutoff[2] - 1) = 1.0; } - tBasisE->orthogonalize_last(-1, true); + //tBasisE->orthogonalize_last(-1, true); } else { @@ -1290,7 +1300,7 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, { (*tBasisE)(i, cutoff[2] - 2) = 1.0; } - tBasisE->orthogonalize_last(cutoff[2] - 1, true); + //tBasisE->orthogonalize_last(cutoff[2] - 1, true); // Add the E offset vector as the last E basis column // and reorthonormalize. @@ -1298,8 +1308,11 @@ void ROM_Sampler::Finalize(Array &cutoff, ROM_Options& input, { (*tBasisE)(i, cutoff[2] - 1) = (*initE)(i); } - tBasisE->orthogonalize_last(-1, true); + //tBasisE->orthogonalize_last(-1, true); } + // Mass induced GS orthonormalization. + MassGramSchmidt(input, 1, tBasisV); + MassGramSchmidt(input, 2, tBasisE); SetupEQP_Force(&Xsnap0, &Vsnap0, &Esnap0, tBasisV, tBasisE, input, sol); @@ -1455,7 +1468,7 @@ ROM_Basis::ROM_Basis(ROM_Options const& input, MPI_Comm comm_, mfH1.SetSize(tH1size); mfL2.SetSize(tL2size); - ReadSolutionBases(input.window, S); + ReadSolutionBases(input, S); if (spaceTime) { @@ -2372,8 +2385,11 @@ int ROM_Basis::SolutionSizeFOM() const return (2*H1size) + L2size; // full size, not true DOF size } -void ROM_Basis::ReadSolutionBases(const int window, BlockVector const* S) +void ROM_Basis::ReadSolutionBases(ROM_Options const& input, + BlockVector const* S) { + int const window = input.window; + // Velocity basis formation. if (hyperreductionSamplingType == eqp_energy) { @@ -2427,9 +2443,11 @@ void ROM_Basis::ReadSolutionBases(const int window, BlockVector const* S) // Add the lifted V solution as the last basis vector. AddLastCol_V(*S); } + // Mass induced GS orthonormalization. + MassGramSchmidt(input, 1, basisV); } - else - { + else + { if (!useVX) { basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + @@ -2484,7 +2502,7 @@ void ROM_Basis::ReadSolutionBases(const int window, BlockVector const* S) { (*basisE)(i, tmp_rdime) = 1.0; } - basisE->orthogonalize_last(-1, true); + //basisE->orthogonalize_last(-1, true); } else { @@ -2515,14 +2533,16 @@ void ROM_Basis::ReadSolutionBases(const int window, BlockVector const* S) { (*basisE)(i, tmp_rdime) = 1.0; } - basisE->orthogonalize_last(rdime-1, true); + //basisE->orthogonalize_last(rdime-1, true); // Add the lifted E solution as the last basis vector. AddLastCol_E(*S); } + // Mass induced GS orthonormalization. + MassGramSchmidt(input, 2, basisE); } - else - { + else + { basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + std::to_string(window) + basisIdentifier, tL2size, rdime); } @@ -2664,7 +2684,8 @@ void ROM_Basis::ReadTemporalBases(const int window) } // f is a full vector, not a true vector -void ROM_Basis::ProjectFOMtoROM(Vector const& f, Vector & r, const bool timeDerivative) +void ROM_Basis::ProjectFOMtoROM(Vector const& f, Vector & r, + const bool timeDerivative) { MFEM_VERIFY(r.Size() == rdimx + rdimv + rdime, ""); MFEM_VERIFY(f.Size() == (2*H1size) + L2size, ""); @@ -2681,25 +2702,69 @@ void ROM_Basis::ProjectFOMtoROM(Vector const& f, Vector & r, const bool timeDeri basisX->transposeMult(*fH1, *rX); - for (int i=0; iGetTrueDofs(mfH1); + for (int i = 0; i < H1size; ++i) + sv[i] = f[H1size + i]; - for (int i=0; iMultMv(sv, Msv); - basisV->transposeMult(*fH1, *rV); + for (int i = 0; i < H1size; ++i) + (*gfH1)(i) = Msv[i]; - for (int i=0; iGetTrueDofs(mfH1); - gfL2->GetTrueDofs(mfL2); + for (int i = 0; i < tH1size; ++i) + (*fH1)(i) = mfH1[i]; - for (int i=0; itransposeMult(*fH1, *rV); + } + else + { + for (int i=0; iGetTrueDofs(mfH1); - basisE->transposeMult(*fL2, *rE); + for (int i=0; itransposeMult(*fH1, *rV); + } + + if (hyperreductionSamplingType == eqp_energy) + { + Vector se(L2size), Mse(L2size); + + for (int i = 0; i < L2size; ++i) + se[i] = f[(2*H1size) + i]; + + lhoper->MultMe(se, Mse); + + for (int i = 0; i < L2size; ++i) + (*gfL2)(i) = Mse[i]; + + gfL2->GetTrueDofs(mfL2); + + for (int i = 0; i < tL2size; ++i) + (*fL2)(i) = mfL2[i]; + + basisE->transposeMult(*fL2, *rE); + } + else + { + for (int i=0; iGetTrueDofs(mfL2); + + for (int i=0; itransposeMult(*fL2, *rE); + } for (int i=0; iorthogonalize_last(-1, true); + //basisV->orthogonalize_last(-1, true); } // f is a full vector, not a true vector @@ -2778,7 +2843,7 @@ void ROM_Basis::AddLastCol_E(BlockVector const& f) for (int i = 0; i < tL2size; i++) (*basisE)(i, rdime-1) = mfL2[i]; - basisE->orthogonalize_last(-1, true); + //basisE->orthogonalize_last(-1, true); } // f is a full vector, not a true vector @@ -3587,38 +3652,42 @@ void ROM_Operator::ComputeReducedMv() // TODO: do I need to enforce MPI rank == 0? else if (hyperreduce && hyperreductionSamplingType == eqp_energy) { - // Form inverse of reduced Mv - invMvROM.SetSize(nv); + // Form inverse of reduced Mv + //invMvROM.SetSize(nv); - const int size_H1 = basis->SolutionSizeH1FOM(); - const int tsize_H1 = H1spaceFOM->GetTrueVSize(); - - Wmat = new CAROM::Matrix(size_H1, nv, true); + const int size_H1 = basis->SolutionSizeH1FOM(); + const int tsize_H1 = H1spaceFOM->GetTrueVSize(); - ParGridFunction gf(H1spaceFOM), gf2(H1spaceFOM); - Vector vj(tsize_H1), vi(tsize_H1); - Vector Mvj(size_H1); + Wmat = new CAROM::Matrix(size_H1, nv, true); - for (int j = 0; j < nv; ++j) - { - basis->GetBasisVectorV(false, j, vj); - gf.SetFromTrueDofs(vj); + ParGridFunction gf(H1spaceFOM); + Vector vj(tsize_H1); - // store jth V basis vector in Wmat(:,j) - for (int i=0; iMultMv(gf, Mvj); + for (int j = 0; j < nv; ++j) + { + basis->GetBasisVectorV(false, j, vj); + gf.SetFromTrueDofs(vj); - for (int i = 0; i < nv; ++i) - { - basis->GetBasisVectorV(false, i, vi); - gf2.SetFromTrueDofs(vi); + // store jth V basis vector in Wmat(:,j) + for (int i = 0; i < size_H1; ++i) + { + (*Wmat)(i,j) = gf[i]; + } - invMvROM(i,j) = gf2 * Mvj; - } - } - invMvROM.Invert(); - } + //operFOM->MultMv(gf, Mvj); + + //for (int i = 0; i < nv; ++i) + //{ + // basis->GetBasisVectorV(false, i, vi); + // gf2.SetFromTrueDofs(vi); + // invMvROM(i,j) = gf2 * Mvj; + //} + } + //invMvROM.Invert(); + } else if (!hyperreduce) { MFEM_ABORT("TODO"); @@ -3668,35 +3737,43 @@ void ROM_Operator::ComputeReducedMe() } // TODO: do I need to enforce MPI rank == 0? else if (hyperreduce && hyperreductionSamplingType == eqp_energy) - { - // Form inverse of reduced Me - invMeROM.SetSize(ne); + { + // Form inverse of reduced Me + //invMeROM.SetSize(ne); - const int size_L2 = basis->SolutionSizeL2FOM(); + const int size_L2 = basis->SolutionSizeL2FOM(); + const int tsize_L2 = L2spaceFOM->GetTrueVSize(); - Wmat_E = new CAROM::Matrix(size_L2, ne, true); + Wmat_E = new CAROM::Matrix(size_L2, ne, true); - Vector ej(size_L2), ei(size_L2); - Vector Mej(size_L2); + ParGridFunction gf(L2spaceFOM); + Vector ej(tsize_L2); - for (int j = 0; j < ne; ++j) - { - basis->GetBasisVectorE(false, j, ej); + //ParGridFunction gf2(L2spaceFOM); + //Vector ei(tsize_L2), Mej(size_L2); - // store jth E basis vector in Wmat_E(:,j) - for (int i=0; iGetBasisVectorE(false, j, ej); + gf.SetFromTrueDofs(ej); - operFOM->MultMe(ej, Mej); + // store jth E basis vector in Wmat_E(:,j) + for (int i = 0; i < size_L2; ++i) + { + (*Wmat_E)(i,j) = gf[i]; + } - for (int i = 0; i < ne; ++i) - { - basis->GetBasisVectorE(false, i, ei); + //operFOM->MultMe(gf, Mej); - invMeROM(i,j) = ei * Mej; - } - } - invMeROM.Invert(); - } + //for (int i = 0; i < ne; ++i) + //{ + // basis->GetBasisVectorE(false, i, ei); + // gf2.SetFromTrueDofs(ei); + // invMeROM(i,j) = gf2 * Mej; + //} + } + //invMeROM.Invert(); + } else if (!hyperreduce) { MFEM_ABORT("TODO"); @@ -3752,7 +3829,9 @@ void ROM_Operator::Mult(const Vector &x, Vector &y) const } } -void ROM_Operator::InducedInnerProduct(const int id1, const int id2, const int var, const int dim, double &ip) +void ROM_Operator::InducedInnerProduct(const int id1, const int id2, + const int var, const int dim, + double &ip) { ip = 0.0; if (use_sample_mesh) @@ -3774,7 +3853,7 @@ void ROM_Operator::InducedInnerProduct(const int id1, const int id2, const int v operSP->MultMe(xj_sp, Mxj_sp); } else - MFEM_ABORT("Invalid input"); + MFEM_ABORT("Invalid variable index."); for (int k=0; kGetTrueVSize(); + int const sizeH1 = input.H1FESpace->GetVSize(); + + MFEM_VERIFY(size == tsizeH1, "Input size mismatch."); + MFEM_VERIFY(basisMat->numRows() == tsizeH1, "Input size mismatch."); + + ParGridFunction gf(input.H1FESpace), gf2(input.H1FESpace); + Vector vj(tsizeH1), vi(tsizeH1); + Vector Mvj(sizeH1); + + for (int i = 0; i < tsizeH1; ++i) + { + vj[i] = (*basisMat)(i, id1); + } + gf.SetFromTrueDofs(vj); + + for (int i = 0; i < tsizeH1; ++i) + { + vi[i] = (*basisMat)(i, id2); + } + gf2.SetFromTrueDofs(vi); + + input.FOMoper->MultMv(gf, Mvj); + inner_prod = gf2 * Mvj; + } + else if (var == 2) // energy + { + int const tsizeL2 = input.L2FESpace->GetTrueVSize(); + int const sizeL2 = input.L2FESpace->GetVSize(); + + MFEM_VERIFY(size == tsizeL2, "Input size mismatch."); + MFEM_VERIFY(basisMat->numRows() == tsizeL2, "Input size mismatch."); + + ParGridFunction gf(input.L2FESpace), gf2(input.L2FESpace); + Vector ej(tsizeL2), ei(tsizeL2); + Vector Mej(sizeL2); + + for (int i = 0; i < tsizeL2; ++i) + { + ej[i] = (*basisMat)(i, id1); + } + gf.SetFromTrueDofs(ej); + + for (int i = 0; i < sizeL2; ++i) + { + ei[i] = (*basisMat)(i, id2); + } + gf2.SetFromTrueDofs(ei); + + input.FOMoper->MultMe(gf, Mej); + inner_prod = gf2 * Mej; + } + else + { + MFEM_ABORT("Invalid variable index."); + } +} + +void MassGramSchmidt(ROM_Options const& input, int const var, + CAROM::Matrix* basisMat) +{ + MFEM_VERIFY(input.hyperreductionSamplingType == eqp_energy, ""); + int const nrows = basisMat->numRows(); + int const ncols = basisMat->numColumns(); + + if (var == 1) // velocity + { + MFEM_VERIFY(nrows == input.H1FESpace->GetTrueVSize(), ""); + } + else if (var == 2) // energy + { + MFEM_VERIFY(nrows == input.L2FESpace->GetVSize(), ""); + } + else + { + MFEM_ABORT("Invalid variable index."); + } + + double factor; + + for (int work = 0; work < ncols; ++work) + { + // Orthogonalize the column (twice currently). + for (int pass = 0; pass < 2; ++pass) + { + for (int col = 0; col < work; ++col) + { + MassInnerProduct(input, var, nrows, basisMat, col, work, + factor); + for (int row = 0; row < nrows; ++row) + { + (*basisMat)(row, work) -= factor * (*basisMat)(row, col); + } + } + } + + // Normalize the column. + MassInnerProduct(input, var, nrows, basisMat, work, work, factor); + + if (factor > 1.0e-15) // zero tolerance + { + double inv_norm = 1.0 / sqrt(factor); + for (int row = 0; row < nrows; ++row) + { + (*basisMat)(row, work) *= inv_norm; + } + } + } +} + +void MassGramSchmidt(ROM_Options const& input, int const var, + CAROM::Matrix* basisMat, CAROM::Matrix* Rmat) +{ + // The Gram-Scmidt algorithm factorizes input marix basisMat + // in the form QR, where + // Q is orthonormal w.r.t. the mass induced inner product, + // R is the upper triangular matrix of factors involved in the + // orthonormalization. + // Input matrix basisMat is overwritten with Q. + + MFEM_VERIFY(input.hyperreductionSamplingType == eqp_energy, ""); + int const nrows = basisMat->numRows(); + int const ncols = basisMat->numColumns(); + + MFEM_VERIFY(Rmat->numRows() == ncols, ""); + MFEM_VERIFY(Rmat->numColumns() == ncols, ""); + + if (var == 1) // velocity + { + MFEM_VERIFY(nrows == input.H1FESpace->GetTrueVSize(), ""); + } + else if (var == 2) // energy + { + MFEM_VERIFY(nrows == input.L2FESpace->GetVSize(), ""); + } + else + { + MFEM_ABORT("Invalid variable index."); + } + + double factor; + + for (int work = 0; work < ncols; ++work) + { + for (int row = 0; row < ncols; ++row) + { + (*Rmat)(row, work) = 0.0; // initialize to zero + } + + // Orthogonalize the column (twice currently). + for (int pass = 0; pass < 2; ++pass) + { + for (int col = 0; col < work; ++col) + { + MassInnerProduct(input, var, nrows, basisMat, col, work, + factor); + (*Rmat)(col, work) += factor; // added in each pass + for (int row = 0; row < nrows; ++row) + { + (*basisMat)(row, work) -= factor * (*basisMat)(row, col); + } + } + } + + // Normalize the column. + MassInnerProduct(input, var, nrows, basisMat, work, work, factor); + + if (factor > 1.0e-15) // zero tolerance + { + (*Rmat)(work, work) = sqrt(factor); + + double inv_norm = 1.0 / (*Rmat)(work, work); + for (int row = 0; row < nrows; ++row) + { + (*basisMat)(row, work) *= inv_norm; + } + } + } +} + diff --git a/rom/laghos_utils.hpp b/rom/laghos_utils.hpp index 8845fbd0..211a1cfa 100644 --- a/rom/laghos_utils.hpp +++ b/rom/laghos_utils.hpp @@ -60,4 +60,16 @@ void readVec(vector &v, std::string file_name); // count the number of lines in a file int countNumLines(std::string file_name); +// Inner product induced by the mass matrix. +void MassInnerProduct(ROM_Options const& input, int const var, int const size, + CAROM::Matrix* basisMat, int const id1, int const id2, + double& inner_prod); + +// Gram-Schmidt orthonormalization w.r.t. the mass inner product. +void MassGramSchmidt(ROM_Options const& input, int const var, + CAROM::Matrix* basisMat); + +void MassGramSchmidt(ROM_Options const& input, int const var, + CAROM::Matrix* basisMat, CAROM::Matrix* Rmat); + #endif // MFEM_LAGHOS_UTILS From 63990ad74f37089144ec028cb7b4d7e6c411b8c9 Mon Sep 17 00:00:00 2001 From: Chris Vales <52826240+cval26@users.noreply.github.com> Date: Tue, 25 Jun 2024 06:13:10 -0700 Subject: [PATCH 38/39] get L directly from libROM's LQ --- rom/laghos_rom.cpp | 50 +++++++++++++++++----------------------------- 1 file changed, 18 insertions(+), 32 deletions(-) diff --git a/rom/laghos_rom.cpp b/rom/laghos_rom.cpp index c9b4e968..08cba534 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -431,59 +431,47 @@ void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, if (nnls_use_lq) { - // Compute Q^T of the LQ factorization of G. - CAROM::Matrix* Qt_ptr; - Qt_ptr = Gt.qr_factorize(); - - CAROM::Matrix Qt(Qt_ptr->getData(), Qt_ptr->numRows(), - Qt_ptr->numColumns(), Qt_ptr->distributed(), - false); - - // Compute L of the factorization; L is lower triangular. - // G = L * Q --> L = G * Q^T - CAROM::Matrix L(Qt.numColumns(), Qt.numColumns(), false); - Gt.transposeMult(Qt, L); - - // Check for nearly linearly dependent Q rows. - // This is achieved by checking the magnitude of the diagonal - // values of L. + // Compute LQ factorization of G. + std::vector> QR; + Gt.qr_factorize(QR); + const CAROM::Matrix & Qt = *QR[0]; // Q transpose + const CAROM::Matrix & R = *QR[1]; // L transpose + + // Check for nearly linearly dependent Q rows, by checking the + // magnitude of the diagonal values of L (equivalently, R). std::vector row_ind; - for (int i = 0; i < L.numRows(); ++i) + for (int i = 0; i < R.numRows(); ++i) { - if (std::abs(L.item(i, i)) < 1e-12) + if (std::abs(R(i, i)) < 1e-12) { row_ind.push_back(i); std::cout << i << ", "; } } - std::cout << "\n"; - std::cout << "Found " << row_ind.size() << " / " << - L.numRows() << " nearly linearly dependent constraints.\n"; + std::cout << "\nFound " << row_ind.size() << " / " << + R.numRows() << " nearly linearly dependent constraints.\n"; // Compute the RHS vector. CAROM::Vector rhs_ub(Qt.numColumns(), false); Qt.transposeMult(w, rhs_ub); // Compute the new RHS tolerance values. - const double delta = 1.0e-11; + double const delta = 1.0e-11; CAROM::Vector delta_new(rhs_ub.dim(), false); for (int i = 0; i < delta_new.dim(); ++i) { - double denominator = (i + 1) * std::abs(L.item(i, i)); + double denominator = (i + 1) * std::abs(R(i, i)); if (std::abs(denominator) < delta) delta_new(i) = 1.0; else delta_new(i) = delta / denominator; + for (int j = i + 1; j < delta_new.dim(); ++j) { - denominator = (j + 1) * std::abs(L.item(j, i)); - - double temp; - if (std::abs(denominator) < delta) - temp = 1.0; - else - temp = delta / denominator; + denominator = (j + 1) * std::abs(R(i, j)); // L(j, i) + double const temp = std::abs(denominator) < delta ? 1.0 : + delta / denominator; if (temp < delta_new(i)) delta_new(i) = temp; @@ -500,8 +488,6 @@ void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, // Call the NNLS solver. nnls.solve_parallel_with_scalapack(Qt, rhs_lb, rhs_ub, sol); - - delete Qt_ptr; } else { From 298da43d1ba29061f36008507efba7fe22174dfc Mon Sep 17 00:00:00 2001 From: Chris Vales <52826240+cval26@users.noreply.github.com> Date: Wed, 24 Jul 2024 09:41:54 -0700 Subject: [PATCH 39/39] include rom makefile changes --- rom/makefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rom/makefile b/rom/makefile index a4ca7b2b..00c9a4e2 100644 --- a/rom/makefile +++ b/rom/makefile @@ -135,10 +135,10 @@ TEST_FILES = tests/basisComparator.cpp tests/fileComparator.cpp tests/solutionCo cd $(