diff --git a/rom/.gitignore b/rom/.gitignore new file mode 100644 index 00000000..ab1ebe02 --- /dev/null +++ b/rom/.gitignore @@ -0,0 +1,8 @@ +dependencies/ +run/ +laghos +merge +*.o +*.out +job*.sh +logfiles/ diff --git a/rom/laghos.cpp b/rom/laghos.cpp index 1aa8c8f9..2abeaa43 100644 --- a/rom/laghos.cpp +++ b/rom/laghos.cpp @@ -342,8 +342,13 @@ 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.AddOption(&romOptions.sampfreq, "-sampfreq", "--samp-freq", + "Snapshot sampling frequency."); + + args.Parse(); if (!args.Good()) { if (mpi.Root()) { @@ -394,8 +399,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.hyperreductionSamplingType == eqp_energy) + romOptions.GramSchmidt = false; romOptions.spaceTimeMethod = getSpaceTimeMethod(spaceTimeMethod); const bool spaceTime = (romOptions.spaceTimeMethod != no_space_time); @@ -997,14 +1009,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; @@ -1168,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) @@ -1228,22 +1244,43 @@ 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. + // 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[romOptions.window] = new ROM_Basis(romOptions, MPI_COMM_WORLD, sFactorX, sFactorV); + basis[0] = 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, 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 { - 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, @@ -1440,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) @@ -1501,7 +1539,9 @@ 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); + + basis[romOptions.window] = new ROM_Basis(romOptions, + MPI_COMM_WORLD, sFactorX, sFactorV, nullptr, S); basis[romOptions.window]->Init(romOptions, *S); if (romOptions.mergeXV) @@ -1573,7 +1613,8 @@ int main(int argc, char *argv[]) { romOper[0]->ApplyHyperreduction(romS); } - double windowEndpoint = 0.0; + + double windowEndpoint = 0.0; double windowOverlapMidpoint = 0.0; for (ti = 1; !last_step; ti++) { @@ -1803,7 +1844,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; @@ -1868,7 +1910,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. @@ -1913,7 +1955,8 @@ int main(int argc, char *argv[]) } else { - sampler->Finalize(cutoff, romOptions); + // Form the reduced bases from the snapshots. + sampler->Finalize(cutoff, romOptions, *S); } if (myid == 0 && romOptions.parameterID == -1) { outfile_twp << windowEndpoint << ", " << cutoff[0] << ", " << cutoff[1] << ", " << cutoff[2]; @@ -1968,6 +2011,24 @@ int main(int argc, char *argv[]) if (rom_online) { + if (romOptions.hyperreductionSamplingType == eqp_energy) + { + timeLoopTimer.Stop(); + 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(); + } + double window_par = t; if (problem == 7) { @@ -1997,11 +2058,33 @@ 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; 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, S); + + 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); @@ -2031,9 +2114,32 @@ int main(int argc, char *argv[]) { basis[romOptions.window]->ProjectFOMtoROM(*S, romS); } - if (myid == 0) + if (myid == 0) + { + cout << "Window " << romOptions.window << ": initial romS norm " << romS.Norml2() << endl; + } + + if (romOptions.hyperreductionSamplingType == eqp_energy) { - 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. + 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]; @@ -2042,7 +2148,7 @@ int main(int argc, char *argv[]) { 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); @@ -2088,6 +2194,24 @@ int main(int argc, char *argv[]) << ",\t|e| = " << setprecision(10) << sqrt(tot_norm) << endl; } + + if (romOptions.hyperreductionSamplingType == eqp_energy) + { + timeLoopTimer.Stop(); + 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(); + } // Make sure all ranks have sent their 'v' solution before initiating // another set of GLVis connections (one from each rank): @@ -2163,6 +2287,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) @@ -2189,9 +2326,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(); @@ -2466,8 +2603,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 7e9a3a9b..08cba534 100644 --- a/rom/laghos_rom.cpp +++ b/rom/laghos_rom.cpp @@ -57,9 +57,9 @@ 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 = generator_X->isNextSample(t); + + Vector dSdt; if (!sns && rhsBasis) { dSdt.SetSize(S.Size()); @@ -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) { @@ -152,7 +152,7 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p } } - const bool sampleE = generator_E->isNextSample(t); + const bool sampleE = generator_E->isNextSample(t); if (sampleE) { @@ -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) { @@ -338,6 +338,81 @@ 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& 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, + 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, ""); + MFEM_VERIFY(v_j_e.Size() == h1dofs_cnt*dim, ""); + } + + if (equationE) + { + MFEM_VERIFY(re.Size() == nqp, ""); + 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 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. + 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 + + trial_fe.CalcShape(ip, eshape); // Energy shape function + + if (equationV) + { + // Inner product, on this element, with the jth V basis vector. + rv[q] = (v_j_e * Vloc_force) * (eshape * unitE); + } + if (equationE) + { + // 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 +} + #include "linalg/NNLS.h" void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, @@ -346,23 +421,88 @@ 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); - CAROM::Vector rhs_lb(rhs_ub); - CAROM::Vector rhs_Gw(rhs_ub); + // TODO: Make it a user input option. + bool const nnls_use_lq = true; - const double delta = 1.0e-11; - for (int i=0; i> 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 < R.numRows(); ++i) + { + if (std::abs(R(i, i)) < 1e-12) + { + row_ind.push_back(i); + std::cout << i << ", "; + } + } + 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. + 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(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(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; + } + } + + // Compute the upper and lower bound RHS vectors. + CAROM::Vector rhs_lb(rhs_ub); + for (int i = 0; i < rhs_ub.dim(); ++i) + { + rhs_lb(i) -= delta_new(i); + rhs_ub(i) += delta_new(i); + } + + // Call the NNLS solver. + nnls.solve_parallel_with_scalapack(Qt, rhs_lb, rhs_ub, sol); } + else + { + CAROM::Vector rhs_ub(rhs_Gw); + CAROM::Vector rhs_lb(rhs_Gw); - //nnls.normalize_constraints(Gt, rhs_lb, rhs_ub); - nnls.solve_parallel_with_scalapack(Gt, rhs_lb, rhs_ub, sol); + const double delta = 1.0e-11; + for (int i=0; iGetWeights(). - - Vector v_i(tH1size); + // 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 + // + // 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); @@ -578,7 +739,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; iGetIntegrationRule(); + 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) x NQ; store its transpose Gt. + // The first NBv * nsnap rows of G hold the velocity constraints; + // the remaining NBe * nsnap rows hold the energy constraints. + CAROM::Matrix Gt(NQ, (NBv + NBe) * nsnap, true); + + // row index of G where energy constraints start + const int estart = NBv * nsnap; + + 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 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; iResetQuadratureData(); + 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 + // + // e_j^T * F(v_i,e_i,x_i) * 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 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 + // + // 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, 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 + // 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 + gf2H1 = S_v; // ith velocity snapshot + + for (int e=0; eGetQuadData(), + v_j_e, 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(), + v_j_e, 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; jnumRows()); + for (int k = 0; k < basisE->numRows(); ++k) + { + Ediff[k] = We(k, j); + } + gfL2.SetFromTrueDofs(Ediff); // jth energy basis vector + gf2H1 = S_v; // ith velocity snapshot + + for (int e=0; eGetQuadData(), + v_j_e, 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(), ""); @@ -608,12 +1031,29 @@ 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(), ""); - SetupEQP_Force_Eq(snapX, snapV, snapE, basisV, basisE, input, false); - SetupEQP_Force_Eq(snapX, snapV, snapE, basisV, basisE, input, true); + 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) + { + // 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, + BlockVector 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(); @@ -713,12 +1153,159 @@ 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; } + 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 + // 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); + + CAROM::Matrix* tBasisV = nullptr; + CAROM::Matrix* tBasisE = nullptr; + + // V basis construction. + if (input.window == 0 && initV->norm() < 1e-15) + { + // 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); + } + + // 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. + for (int i = 0; i < tL2size; ++i) + { + (*tBasisE)(i, cutoff[2] - 1) = 1.0; + } + //tBasisE->orthogonalize_last(-1, true); + } + else + { + 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) + { + (*tBasisE)(i, j) = (*basisE)(i, j); + } + } + // Add the energy identity as the penultimate E basis column + // and reorthonormalize. + for (int i = 0; i < tL2size; ++i) + { + (*tBasisE)(i, cutoff[2] - 2) = 1.0; + } + //tBasisE->orthogonalize_last(cutoff[2] - 1, true); + + // Add the E offset vector as the last E basis column + // and reorthonormalize. + for (int i = 0; i < tL2size; ++i) + { + (*tBasisE)(i, cutoff[2] - 1) = (*initE)(i); + } + //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); + + delete tBasisV, tBasisE; + } + delete generator_X; delete generator_V; delete generator_E; @@ -738,11 +1325,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, ""); @@ -805,8 +1392,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), @@ -865,7 +1454,8 @@ ROM_Basis::ROM_Basis(ROM_Options const& input, MPI_Comm comm_, const double sFac mfH1.SetSize(tH1size); mfL2.SetSize(tL2size); - ReadSolutionBases(input.window); + ReadSolutionBases(input, S); + if (spaceTime) { ReadTemporalBases(input.window); @@ -1712,6 +2302,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 @@ -1733,6 +2326,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); @@ -1776,17 +2371,208 @@ 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(ROM_Options const& input, + BlockVector const* S) { - if (!useVX) - basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + std::to_string(window) + basisIdentifier, tH1size, rdimv); + int const window = input.window; + + // Velocity basis formation. + if (hyperreductionSamplingType == eqp_energy) + { + 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"); + + 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); + } + // Mass induced GS orthonormalization. + MassGramSchmidt(input, 1, basisV); + } + else + { + if (!useVX) + { + basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + + std::to_string(window) + basisIdentifier, tH1size, rdimv); + } + } + + // Energy basis formation. + 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; - basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + std::to_string(window) + basisIdentifier, tL2size, rdime); + // 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); - if (useXV) - basisX = basisV; + // Add the lifted E solution as the last basis vector. + AddLastCol_E(*S); + } + // Mass induced GS orthonormalization. + MassGramSchmidt(input, 2, basisE); + } + 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); + + 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; @@ -1808,7 +2594,7 @@ void ROM_Basis::ReadSolutionBases(const int window) 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; @@ -1822,7 +2608,7 @@ void ROM_Basis::ReadSolutionBases(const int window) 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; @@ -1845,7 +2631,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 { @@ -1884,7 +2670,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, ""); @@ -1901,25 +2688,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; itransposeMult(*fL2, *rE); + gfH1->GetTrueDofs(mfH1); + + 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; iGetTrueDofs(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, ""); + + for (int i=0; iGetTrueDofs(mfH1); + + for (int i = 0; i < tH1size; i++) + (*basisV)(i, rdimv-1) = mfH1[i]; + + //basisV->orthogonalize_last(-1, true); +} + +// f is a full vector, not a true vector +void ROM_Basis::AddLastCol_E(BlockVector 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_last(-1, true); +} + // f is a full vector, not a true vector void ROM_Basis::LiftROMtoFOM(Vector const& r, Vector & f) { @@ -2353,22 +3232,41 @@ 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(); ComputeReducedMe(); } - if (hyperreduce && hyperreductionSamplingType == eqp) - { - // TODO: are reduced mass matrices needed for EQP, or just W matrices? - - ComputeReducedMv(); - ComputeReducedMe(); - - ReadSolutionNNLS(input, "run/nnlsV", eqpI, eqpW); - ReadSolutionNNLS(input, "run/nnlsE", eqpI_E, eqpW_E); - } + if (hyperreduce && hyperreductionSamplingType == eqp) + { + // basic EQP + + // 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); + } + else if (hyperreduce && hyperreductionSamplingType == eqp_energy) + { + // energy-conserving EQP + + // Computes Phi_v^T * Mv * Phi_v (similarly for energy) needed to + // form the RHS vectors to time march. + 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); + } } void ROM_Operator::ReadSolutionNNLS(ROM_Options const& input, string basename, @@ -2736,6 +3634,45 @@ void ROM_Operator::ComputeReducedMv() for (int i=0; iSolutionSizeH1FOM(); + const int tsize_H1 = H1spaceFOM->GetTrueVSize(); + + Wmat = new CAROM::Matrix(size_H1, nv, true); + + ParGridFunction gf(H1spaceFOM); + Vector vj(tsize_H1); + + //ParGridFunction gf2(H1spaceFOM); + //Vector vi(tsize_H1), Mvj(size_H1); + + for (int j = 0; j < nv; ++j) + { + basis->GetBasisVectorV(false, j, vj); + gf.SetFromTrueDofs(vj); + + // store jth V basis vector in Wmat(:,j) + for (int i = 0; i < size_H1; ++i) + { + (*Wmat)(i,j) = gf[i]; + } + + //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) { @@ -2783,6 +3720,45 @@ void ROM_Operator::ComputeReducedMe() for (int i=0; iSolutionSizeL2FOM(); + const int tsize_L2 = L2spaceFOM->GetTrueVSize(); + + Wmat_E = new CAROM::Matrix(size_L2, ne, true); + + ParGridFunction gf(L2spaceFOM); + Vector ej(tsize_L2); + + //ParGridFunction gf2(L2spaceFOM); + //Vector ei(tsize_L2), Mej(size_L2); + + for (int j = 0; j < ne; ++j) + { + basis->GetBasisVectorE(false, j, ej); + gf.SetFromTrueDofs(ej); + + // store jth E basis vector in Wmat_E(:,j) + for (int i = 0; i < size_L2; ++i) + { + (*Wmat_E)(i,j) = gf[i]; + } + + //operFOM->MultMe(gf, Mej); + + //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) { @@ -2839,7 +3815,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) @@ -2861,7 +3839,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; kSetPointsEQP(eqpI_E); } -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(); @@ -4027,14 +5006,10 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res) const const int nqe = ir->GetWeights().Size(); - DenseMatrix vshape, loc_force; - Vector shape, unitE, rhs; - + DenseMatrix grad_vshape, loc_force; Array vdofs; - - Vector v_e; - - res = 0.0; + + res = 0.0; int eprev = -1; int dof = 0; @@ -4126,17 +5101,15 @@ 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); + grad_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); + // 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. @@ -4145,40 +5118,64 @@ void ROM_Operator::ForceIntegratorEQP(Vector & res) const 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); } } } - loc_force *= eqpW[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()); - - unitE.SetSize(shape.Size()); - unitE = 1.0; - - rhs.SetSize(h1dofs_cnt * dim); - - v_e.SetSize(rhs.Size()); - - const int eos = elemIndex * nvdof; - for (int j=0; jCalcShape(ip, eshape); + unitE = 1.0; + + const int eos = elemIndex * nvdof; + + if (energy_conserve) + { + // energy-conserving EQP + for (int j = 0; j < rdim; ++j) + { + // 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 jth V basis vector. + res[j] += (v_e * Vloc_force) * (eshape * unitE); + } + } + else + { + // basic EQP + 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 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 + // 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); } -void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res) const +void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res, + bool energy_conserve) const { const IntegrationRule *ir = operFOM->GetIntegrationRule(); const int rdim = basis->GetDimE(); @@ -4187,12 +5184,10 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res) const const int nqe = ir->GetWeights().Size(); - DenseMatrix vshape, loc_force; - Vector shape, rhs; - + DenseMatrix grad_vshape, loc_force; Array vdofs, edofs; - Vector v_e, w_e; + Vector v_e; res = 0.0; @@ -4301,13 +5296,11 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res) const const int l2dofs_cnt = trial_fe->GetDof(); - vshape.SetSize(h1dofs_cnt, dim); + grad_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); + // 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. @@ -4316,49 +5309,76 @@ void ROM_Operator::ForceIntegratorEQP_E(Vector const& v, Vector & res) const 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); } } } - 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()); - - rhs.SetSize(h1dofs_cnt * dim); - - w_e.SetSize(nedof); - - const int eos = elemIndex * nedof; - for (int j=0; jCalcShape(ip, eshape); + + const int eos = elemIndex * nedof; + + if (energy_conserve) + { + // energy-conserving EQP + for (int j = 0; j < rdim; j++) + { + // w_e: jth E basis vector's DOFs on this element + for (int k=0; kGetDimV()); - 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()); - ForceIntegratorEQP_E(v, res); + ForceIntegratorEQP_E(v, res, energy_conserve); basis->LiftROMtoFOM_dEdt(res, rhs); } @@ -4375,10 +5395,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(); @@ -4398,36 +5417,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) @@ -4443,3 +5467,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 dd20fa1e..8d63194e 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"); @@ -257,7 +259,11 @@ 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 + + // snapshot sampling frequency (sample every sampfreq timestep) + int sampfreq = 1; }; static double* getGreedyParam(ROM_Options& romOptions, const char* greedyParam) @@ -479,13 +485,14 @@ 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), + hyperreductionSamplingType(input.hyperreductionSamplingType) { const int window = input.window; // 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; @@ -625,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, BlockVector const& sol); int MaxNumSamples() { @@ -683,6 +690,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. @@ -817,11 +826,15 @@ 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, + 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, 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 @@ -831,13 +844,15 @@ 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(ROM_Options const& input, + BlockVector const* S = nullptr); void ReadTemporalBases(const int window); void ProjectFOMtoROM(Vector const& f, Vector & r, @@ -846,6 +861,10 @@ class ROM_Basis void ProjectFOMtoROM_V(Vector const& f, Vector & r, const bool timeDerivative=false); + 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); void LiftROMtoFOM_dEdt(Vector const& r, Vector & f); @@ -976,6 +995,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; @@ -1222,11 +1246,15 @@ 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, + bool energy_conserve = false) const; + void ForceIntegratorEQP_E(Vector const& v, 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; + HyperreductionSamplingType getSamplingType() const; void InitEQP() const; diff --git a/rom/laghos_solver.cpp b/rom/laghos_solver.cpp index 4b04f73b..368e42a0 100644 --- a/rom/laghos_solver.cpp +++ b/rom/laghos_solver.cpp @@ -279,7 +279,7 @@ void LagrangianHydroOperator::SolveVelocity(const Vector &S, if (use_eqp) { - rom_op->ForceIntegratorEQP_FOM(rhs); + rom_op->ForceIntegratorEQP_FOM(rhs, rom_op->getSamplingType() == eqp_energy); } else { @@ -410,7 +410,7 @@ void LagrangianHydroOperator::SolveEnergy(const Vector &S, const Vector &v, timer.sw_force.Start(); if (use_eqp) { - rom_op->ForceIntegratorEQP_E_FOM(v, e_rhs); + rom_op->ForceIntegratorEQP_E_FOM(v, e_rhs, rom_op->getSamplingType() == eqp_energy); } else { diff --git a/rom/laghos_utils.cpp b/rom/laghos_utils.cpp index fc5b14e4..f06a9422 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) @@ -507,11 +510,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; @@ -865,3 +876,193 @@ void PrintL2NormsOfParGridFunctions(const int rank, const std::string& name, Par cout << rank << ": " << name << " DIFF norm " << sqrt(diffglob2) << endl; cout << rank << ": " << name << " Rel. DIFF norm " << sqrt(diffglob2)/sqrt(fomglob2) << endl; } + +void MassInnerProduct(ROM_Options const& input, int const var, int const size, + CAROM::Matrix* basisMat, int const id1, int const id2, + double& inner_prod) +{ + MFEM_VERIFY(input.hyperreductionSamplingType == eqp_energy, ""); + inner_prod = 0.0; + + if (var == 1) // velocity + { + int const tsizeH1 = input.H1FESpace->GetTrueVSize(); + 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 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 $(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(); }