diff --git a/examples/prom/mixed_nonlinear_diffusion_eqp.cpp b/examples/prom/mixed_nonlinear_diffusion_eqp.cpp new file mode 100644 index 000000000..6fc797117 --- /dev/null +++ b/examples/prom/mixed_nonlinear_diffusion_eqp.cpp @@ -0,0 +1,1006 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// Functions used by mixed_nonlinear_diffusion.cpp with EQP. + +#include "mfem/Utilities.hpp" + +using namespace mfem; +using namespace std; + +///#include "linalg/NNLS.h" +//#include "mfem/SampleMesh.hpp" +#include "mixed_nonlinear_diffusion_eqp.hpp" + +// static bool nonlinear_problem; +// static double diffusion_c; + +/*double NonlinearCoefficient(const double p); + +double NonlinearCoefficient(const double p) +{ + if (nonlinear_problem) + return diffusion_c + p; + else + return 1.0; +} +*/ + +void AssembleElementMatrix_VectorFEMassIntegrator(Coefficient *Q, + const FiniteElement &el, + ElementTransformation &Trans, + DenseMatrix &elmat) +{ + int dof = el.GetDof(); + int spaceDim = Trans.GetSpaceDim(); + + double w; + + DenseMatrix trial_vshape(dof, spaceDim); + + elmat.SetSize(dof); + elmat = 0.0; + + const IntegrationRule *ir = NULL; + if (ir == NULL) + { + // int order = 2 * el.GetOrder(); + int order = Trans.OrderW() + 2 * el.GetOrder(); + ir = &IntRules.Get(el.GetGeomType(), order); + } + + for (int i = 0; i < ir->GetNPoints(); i++) + { + const IntegrationPoint &ip = ir->IntPoint(i); + + Trans.SetIntPoint(&ip); + + el.CalcVShape(Trans, trial_vshape); + + w = ip.weight * Trans.Weight(); + + if (Q) + { + w *= Q -> Eval (Trans, ip); + } + AddMult_a_AAt (w, trial_vshape, elmat); + } +} + + +void GetEQPCoefficients_VectorFEMassIntegrator(ParFiniteElementSpace *fesR, + std::vector const& rw, std::vector const& qp, + const IntegrationRule *ir, + CAROM::Matrix const& V, Vector & res) +{ + const int rdim = V.numColumns(); + MFEM_VERIFY(V.numRows() == fesR->GetTrueVSize(), ""); + MFEM_VERIFY(rw.size() == qp.size(), ""); + + const int ne = fesR->GetNE(); + const int nqe = ir->GetWeights().Size(); + + ElementTransformation *eltrans; + DofTransformation * doftrans; + const FiniteElement *fe = NULL; + Array vdofs; + + DenseMatrix trial_vshape; + + Vector Vx; + res.SetSize(rdim * rdim * rw.size()); + res = 0.0; + + // For the parallel case, we must get all DOFs of V on sampled elements. + CAROM::Matrix Vs; + + // Since V only has rows for true DOFs, we use a ParGridFunction to get all DOFs. + int eprev = -1; + int dof = 0; + int elemCount = 0; + + // First, find all sampled elements. + for (int i=0; iGetElementVDofs(e, vdofs); + if (dof > 0) + { + MFEM_VERIFY(dof == vdofs.Size(), "All elements must have same DOF size"); + } + dof = vdofs.Size(); + eprev = e; + elemCount++; + } + } + + // Now set Vs. + // Note that, ideally, these FOM data structures and operations should be + // avoided, but this only done in setup. + Vs.setSize(elemCount * dof, rdim); + ParGridFunction v_gf(fesR); + Vector vtrue(fesR->GetTrueVSize()); + + for (int j=0; jGetElementVDofs(e, vdofs); + + for (int k=0; k= 0) ? vdofs[k] : -1 - vdofs[k]; + const double vk = (vdofs[k] >= 0) ? v_gf[dofk] : -v_gf[dofk]; + Vs((elemCount * dof) + k, j) = vk; + } + + eprev = e; + elemCount++; + } + } + } + + eprev = -1; + elemCount = 0; + int spaceDim = 0; + + for (int i=0; iIntPoint(qpi); + + if (e != eprev) // Update element transformation + { + doftrans = fesR->GetElementVDofs(e, vdofs); + fe = fesR->GetFE(e); + eltrans = fesR->GetElementTransformation(e); + + if (doftrans) + { + MFEM_ABORT("TODO"); + } + + dof = fe->GetDof(); + spaceDim = eltrans->GetSpaceDim(); + trial_vshape.SetSize(dof, spaceDim); + Vx.SetSize(spaceDim); + + eprev = e; + elemCount++; + } + + // Integrate at the current point + + eltrans->SetIntPoint(&ip); + fe->CalcVShape(*eltrans, trial_vshape); + + double w = eltrans->Weight() * rw[i]; // using rw[i] instead of ip.weight + + // Note that the coefficient Q is omitted: w *= Q -> Eval(*eltrans, ip); + + for (int jx=0; jx const& rw, std::vector const& qp, + const IntegrationRule *ir, Coefficient *Q, + CAROM::Matrix const& V, CAROM::Vector const& x, const int rank, Vector & res) +{ + const int rdim = V.numColumns(); + MFEM_VERIFY(rw.size() == qp.size(), ""); + MFEM_VERIFY(x.dim() == rdim, ""); + MFEM_VERIFY(V.numRows() == fesR->GetTrueVSize(), ""); + + MFEM_VERIFY(rank == 0, + "TODO: generalize to parallel. This uses full dofs in V, which has true dofs"); + + const int nqe = ir->GetWeights().Size(); + + ElementTransformation *eltrans; + DofTransformation * doftrans; + const FiniteElement *fe = NULL; + Array vdofs; + + DenseMatrix trial_vshape; + + Vector Vx; + res.SetSize(rdim); + res = 0.0; + + int eprev = -1; + int dof = 0; + int spaceDim = 0; + + // Note that the alternative version VectorFEMassIntegrator_ComputeReducedEQP_Fast + // of this function is optimized by storing some intermediate computations. + + for (int i=0; iIntPoint(qpi); + + if (e != eprev) // Update element transformation + { + doftrans = fesR->GetElementVDofs(e, vdofs); + fe = fesR->GetFE(e); + eltrans = fesR->GetElementTransformation(e); + + if (doftrans) + { + MFEM_ABORT("TODO"); + } + + dof = fe->GetDof(); + spaceDim = eltrans->GetSpaceDim(); + trial_vshape.SetSize(dof, spaceDim); + Vx.SetSize(spaceDim); + + eprev = e; + } + + // Integrate at the current point + + eltrans->SetIntPoint(&ip); + fe->CalcVShape(*eltrans, trial_vshape); + + double w = eltrans->Weight() * rw[i]; // using rw[i] instead of ip.weight + + if (Q) + { + w *= Q -> Eval(*eltrans, ip); + } + + // Lift Vx at ip. + Vx = 0.0; + for (int k=0; k= 0) ? vdofs[k] : -1 - vdofs[k]; + + double Vx_k = 0.0; + for (int j=0; j= 0) ? vdofs[l] : -1 - vdofs[l]; + const double s = (vdofs[l] >= 0) ? 1.0 : -1.0; + Vjk += s * V(dofl, j) * trial_vshape(l, k); + } + + rj += Vx[k] * Vjk; + } + + res[j] += w * rj; + } + } +} + +void VectorFEMassIntegrator_ComputeReducedEQP_Fast(ParFiniteElementSpace *fesR, + std::vector const& qp, const IntegrationRule *ir, + Coefficient *Q, CAROM::Vector const& x, + Vector const& coef, Vector & res) +{ + const int rdim = x.dim(); + const int nqe = ir->GetWeights().Size(); + ElementTransformation *eltrans; + + res.SetSize(rdim); + res = 0.0; + + int eprev = -1; + + for (int i=0; iIntPoint(qpi); + + if (e != eprev) // Update element transformation + { + eltrans = fesR->GetElementTransformation(e); + eprev = e; + } + + eltrans->SetIntPoint(&ip); + const double q_ip = Q -> Eval(*eltrans, ip); + + for (int j=0; j const& rw, std::vector const& qp, + const IntegrationRule *ir, CAROM::Matrix const& V, + Vector & res) +{ + // Assumption: fes is an L2 space, where true DOFs and full DOFs are the same. + // This allows for using V(dof,j), where dof is from fes->GetElementVDofs. + MFEM_VERIFY(fes->GetTrueVSize() == fes->GetTrueVSize(), ""); + + const int rdim = V.numColumns(); + MFEM_VERIFY(rw.size() == qp.size(), ""); + MFEM_VERIFY(V.numRows() == fes->GetTrueVSize(), ""); + + const int nqe = ir->GetWeights().Size(); + + ElementTransformation *eltrans; + DofTransformation * doftrans; + const FiniteElement *fe = NULL; + Array vdofs; + + Vector shape; + + int eprev = -1; + int dof = 0; + + for (int i=0; iIntPoint(qpi); + + if (e != eprev) // Update element transformation + { + doftrans = fes->GetElementVDofs(e, vdofs); + fe = fes->GetFE(e); + eltrans = fes->GetElementTransformation(e); + + if (doftrans) + { + MFEM_ABORT("TODO"); + } + + dof = fe->GetDof(); + + shape.SetSize(dof); + + eprev = e; + + if (i == 0) + { + res.SetSize(rdim * rw.size() * dof); + res = 0.0; + } + } + + // Integrate at the current point + + eltrans->SetIntPoint(&ip); + + fe->CalcPhysShape(*eltrans, shape); + + double w = eltrans->Weight() * rw[i]; // using rw[i] instead of ip.weight + + for (int j=0; j= 0) ? vdofs[l] : -1 - vdofs[l]; + const double s = (vdofs[l] >= 0) ? 1.0 : -1.0; + res[l + (j * dof) + (i * rdim * dof)] = w * s * V(dofl, + j) * shape(l) * shape(l); + } + } + } +} + +void LinearMassIntegrator_ComputeReducedEQP(ParFiniteElementSpace *fes, + std::vector const& rw, std::vector const& qp, + const IntegrationRule *ir, Coefficient *Q, + CAROM::Matrix const& V, CAROM::Vector & res) +{ + const int rdim = V.numColumns(); + MFEM_VERIFY(rw.size() == qp.size(), ""); + MFEM_VERIFY(res.dim() == rdim, ""); + MFEM_VERIFY(V.numRows() == fes->GetTrueVSize(), ""); + + ParGridFunction f_gf(fes); + f_gf.ProjectCoefficient(*Q); + Vector sv(fes->GetTrueVSize()); + f_gf.GetTrueDofs(sv); + + const int nqe = ir->GetWeights().Size(); + + ElementTransformation *eltrans; + DofTransformation * doftrans; + const FiniteElement *fe = NULL; + Array vdofs; + + Vector shape; + + res = 0.0; + + int eprev = -1; + int dof = 0; + + // Note that the alternative version LinearMassIntegrator_ComputeReducedEQP_Fast + // of this function is optimized by storing some intermediate computations. + + for (int i=0; iIntPoint(qpi); + + if (e != eprev) // Update element transformation + { + doftrans = fes->GetElementVDofs(e, vdofs); + fe = fes->GetFE(e); + eltrans = fes->GetElementTransformation(e); + + if (doftrans) + { + MFEM_ABORT("TODO"); + } + + dof = fe->GetDof(); + + shape.SetSize(dof); + + eprev = e; + } + + // Integrate at the current point + + eltrans->SetIntPoint(&ip); + + fe->CalcPhysShape(*eltrans, shape); + + double w = eltrans->Weight() * rw[i]; // using rw[i] instead of ip.weight + + //w *= Q -> Eval(*eltrans, ip); + + for (int j=0; j= 0) ? vdofs[l] : -1 - vdofs[l]; + const double s = (vdofs[l] >= 0) ? 1.0 : -1.0; + //Vj += s * V(dofl, j) * shape(l); + Vj += s * V(dofl, j) * shape(l) * sv[dofl] * shape(l); + } + + res(j) += w * Vj; + } + } +} + +void LinearMassIntegrator_ComputeReducedEQP_Fast(ParFiniteElementSpace *fes, + ParGridFunction & f_gf, + std::vector const& rw, std::vector const& qp, + const IntegrationRule *ir, Coefficient *Q, + CAROM::Matrix const& V, Vector const& coef, + CAROM::Vector & res) +{ + const int rdim = V.numColumns(); + MFEM_VERIFY(rw.size() == qp.size(), ""); + MFEM_VERIFY(V.numRows() == fes->GetTrueVSize(), ""); + + const int nqe = ir->GetWeights().Size(); + + DofTransformation * doftrans; + const FiniteElement *fe = NULL; + Array vdofs; + + res = 0.0; + + int eprev = -1; + int dof = 0; + + for (int i=0; iIntPoint(qpi); + + if (e != eprev) // Update element transformation + { + doftrans = fes->GetElementVDofs(e, vdofs); + fe = fes->GetFE(e); + + if (doftrans) + { + MFEM_ABORT("TODO"); + } + + dof = fe->GetDof(); + + eprev = e; + + f_gf.ProjectCoefficient(*Q, vdofs); + + if (i == 0) + { + MFEM_VERIFY(coef.Size() == rdim * rw.size() * dof, ""); + } + } + + // Integrate at the current point + + for (int j=0; j= 0) ? vdofs[l] : -1 - vdofs[l]; + Vj += f_gf[dofl] * coef[l + (j * dof) + (i * rdim * dof)]; + } + + res(j) += Vj; + } + } +} + +void ComputeElementRowOfG(const IntegrationRule *ir, Array const& vdofs, + Coefficient *Q, Vector const& u, Vector const& v, + FiniteElement const& fe, ElementTransformation & Trans, Vector & r) +{ + MFEM_VERIFY(r.Size() == ir->GetNPoints(), ""); + int dof = fe.GetDof(); + int spaceDim = Trans.GetSpaceDim(); + + Vector u_i(spaceDim); + Vector v_i(spaceDim); + + DenseMatrix trial_vshape(dof, spaceDim); + + for (int i = 0; i < ir->GetNPoints(); i++) + { + const IntegrationPoint &ip = ir->IntPoint(i); + + Trans.SetIntPoint(&ip); + + fe.CalcVShape(Trans, trial_vshape); + + double w = Trans.Weight(); + + u_i = 0.0; + v_i = 0.0; + + for (int j=0; j= 0) ? vdofs[j] : -1 - vdofs[j]; + const double s = (vdofs[j] >= 0) ? 1.0 : -1.0; + for (int k=0; k Eval (Trans, ip); + } + + r[i] = 0.0; + for (int k=0; k const& vdofs, + Vector const& s, Vector const& p, + FiniteElement const& fe, ElementTransformation & Trans, Vector & r) +{ + MFEM_VERIFY(r.Size() == ir->GetNPoints(), ""); + int dof = fe.GetDof(); + + Vector shape(dof); + + for (int i = 0; i < ir->GetNPoints(); i++) + { + const IntegrationPoint &ip = ir->IntPoint(i); + + Trans.SetIntPoint(&ip); + + fe.CalcPhysShape(Trans, shape); + + double w = Trans.Weight(); + + double u_i = 0.0; + double v_i = 0.0; + + for (int j=0; j= 0) ? vdofs[j] : -1 - vdofs[j]; + const double sd = (vdofs[j] >= 0) ? 1.0 : -1.0; + u_i += sd * p[dofj] * shape(j); + v_i += sd * s[dofj] * shape(j); + } + + r[i] = w * u_i * v_i; + } +} + +void PreconditionNNLS(ParFiniteElementSpace *fespace, + BilinearFormIntegrator *massInteg, + const CAROM::Matrix* B, + const int snapshot, + CAROM::Matrix & Gt) +{ + const int NB = B->numColumns(); + const int NQ = Gt.numRows(); + const int nsnap = Gt.numColumns() / NB; + + ParBilinearForm M(fespace); + + M.AddDomainIntegrator(massInteg); + M.Assemble(); + M.Finalize(); + HypreParMatrix *Mmat = M.ParallelAssemble(); + + CAROM::Matrix Mhat(NB, NB, false); + ComputeCtAB(*Mmat, *B, *B, Mhat); + Mhat.inverse(); + + CAROM::Vector PG(NB, false); + + for (int i = (snapshot >= 0 ? snapshot : 0); + i < (snapshot >= 0 ? snapshot+1 : nsnap); ++i) + { + for (int m=0; mGetNPoints(); + const int ne = fespace_R->GetNE(); + const int NB = BR->numColumns(); + const int NQ = ne * nqe; + const int nsnap = BR_snapshots->numColumns(); + + MFEM_VERIFY(nsnap == BW_snapshots->numColumns() || + nsnap + nsets == BW_snapshots->numColumns(), ""); + MFEM_VERIFY(BR->numRows() == BR_snapshots->numRows(), ""); + MFEM_VERIFY(BR->numRows() == fespace_R->GetTrueVSize(), ""); + MFEM_VERIFY(BW_snapshots->numRows() == fespace_W->GetTrueVSize(), ""); + + const bool skipFirstW = (nsnap + nsets == BW_snapshots->numColumns()); + + // Compute G of size (NB * nsnap) x NQ, but only store its transpose Gt. + CAROM::Matrix Gt(NQ, NB * nsnap, true); + + // For 0 <= j < NB, 0 <= i < nsnap, 0 <= e < ne, 0 <= m < nqe, + // G(j + (i*NB), (e*nqe) + m) + // is the coefficient of v_j^T M(p_i) V 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 p_i(BW_snapshots->numRows()); + Vector v_i(BR_snapshots->numRows()); + Vector v_j(BR->numRows()); + + Vector r(nqe); + + int skip = 0; + const int nsnapPerSet = nsnap / nsets; + if (skipFirstW) + { + MFEM_VERIFY(nsets * nsnapPerSet == nsnap, ""); + skip = 1; + } + + for (int i=0; inumRows(); ++j) + p_i[j] = (*BW_snapshots)(j, i + skip); + + for (int j = 0; j < BR_snapshots->numRows(); ++j) + v_i[j] = (*BR_snapshots)(j, i); + + if (skipFirstW && i > 0 && i % nsnapPerSet == 0) + skip++; + + // Set grid function for a(p) + ParGridFunction p_gf(fespace_W); + + p_gf.SetFromTrueDofs(p_i); + + GridFunctionCoefficient p_coeff(&p_gf); + TransformedCoefficient a_coeff(&p_coeff, NonlinearCoefficient); + + ParGridFunction vi_gf(fespace_R); + vi_gf.SetFromTrueDofs(v_i); + + for (int j=0; jnumRows(); ++k) + v_j[k] = (*BR)(k, j); + + ParGridFunction vj_gf(fespace_R); + vj_gf.SetFromTrueDofs(v_j); + + // TODO: is it better to make the element loop the outer loop? + for (int e=0; e vdofs; + DofTransformation *doftrans = fespace_R->GetElementVDofs(e, vdofs); + const FiniteElement &fe = *fespace_R->GetFE(e); + ElementTransformation *eltrans = fespace_R->GetElementTransformation(e); + + ComputeElementRowOfG(ir0, vdofs, &a_coeff, vi_gf, vj_gf, fe, *eltrans, r); + + for (int m=0; m const& w_el = ir0->GetWeights(); + MFEM_VERIFY(w_el.Size() == nqe, ""); + + CAROM::Vector w(ne * nqe, true); + + for (int i=0; iGetNPoints(); + const int ne = fespace_W->GetNE(); + const int NB = BW->numColumns(); + const int NQ = ne * nqe; + const int nsnap = BS_snapshots->numColumns(); + const int nrows = BW->numRows(); + + MFEM_VERIFY(BS_snapshots->numRows() == BW->numRows(), ""); + + // Compute G of size (NB * nsnap) x NQ, but only store its transpose Gt. + CAROM::Matrix Gt(NQ, NB * nsnap, true); + + Vector s_i(nrows); + Vector p_j(nrows); + + Vector r(nqe); + + for (int i=0; inumRows(); ++j) + s_i[j] = (*BS_snapshots)(j, i); + + for (int j=0; jnumRows(); ++k) + p_j[k] = (*BW)(k, j); + + // TODO: is it better to make the element loop the outer loop? + for (int e=0; e vdofs; + DofTransformation *doftrans = fespace_W->GetElementVDofs(e, vdofs); + const FiniteElement &fe = *fespace_W->GetFE(e); + ElementTransformation *eltrans = fespace_W->GetElementTransformation(e); + + ComputeElementRowOfG_Source(ir0, vdofs, s_i, p_j, fe, *eltrans, r); + + for (int m=0; m const& w_el = ir0->GetWeights(); + MFEM_VERIFY(w_el.Size() == nqe, ""); + + CAROM::Vector w(ne * nqe, true); + + for (int i=0; i elements; + Vector elemCount(pmesh->GetNE()); + elemCount = 0.0; + + for (int i=0; i 1.0e-12) + { + const int e = i / nqe; // Element index + elements.insert(e); + elemCount[e] += 1.0; + } + } + + // Empty sets, since EQP on samples inside elements. + std::set faces, edges, vertices; + CAROM::SampleVisualization(pmesh, elements, elements, faces, edges, + vertices, "EQPvis", &elemCount); +} + diff --git a/examples/prom/mixed_nonlinear_diffusion_eqp.hpp b/examples/prom/mixed_nonlinear_diffusion_eqp.hpp index 0a781f5e0..992b65fa2 100644 --- a/examples/prom/mixed_nonlinear_diffusion_eqp.hpp +++ b/examples/prom/mixed_nonlinear_diffusion_eqp.hpp @@ -15,777 +15,67 @@ using namespace mfem; using namespace std; -// Element matrix assembly, copying element loop from BilinearForm::Assemble(int skip_zeros) -// and element matrix computation from VectorFEMassIntegrator::AssembleElementMatrix. -void AssembleElementMatrix_VectorFEMassIntegrator(Coefficient *Q, - const FiniteElement &el, - ElementTransformation &Trans, - DenseMatrix &elmat) -{ - int dof = el.GetDof(); - int spaceDim = Trans.GetSpaceDim(); - - double w; - - DenseMatrix trial_vshape(dof, spaceDim); - - elmat.SetSize(dof); - elmat = 0.0; - - const IntegrationRule *ir = NULL; - if (ir == NULL) - { - // int order = 2 * el.GetOrder(); - int order = Trans.OrderW() + 2 * el.GetOrder(); - ir = &IntRules.Get(el.GetGeomType(), order); - } - - for (int i = 0; i < ir->GetNPoints(); i++) - { - const IntegrationPoint &ip = ir->IntPoint(i); - - Trans.SetIntPoint(&ip); - - el.CalcVShape(Trans, trial_vshape); - - w = ip.weight * Trans.Weight(); - - if (Q) - { - w *= Q -> Eval (Trans, ip); - } - AddMult_a_AAt (w, trial_vshape, elmat); - } -} +#ifndef MND_EQP +#define MND_EQP // Compute coefficients of the reduced integrator with respect to inputs Q and x // in VectorFEMassIntegrator_ComputeReducedEQP. void GetEQPCoefficients_VectorFEMassIntegrator(ParFiniteElementSpace *fesR, std::vector const& rw, std::vector const& qp, const IntegrationRule *ir, - CAROM::Matrix const& V, Vector & res) -{ - const int rdim = V.numColumns(); - MFEM_VERIFY(V.numRows() == fesR->GetTrueVSize(), ""); - MFEM_VERIFY(rw.size() == qp.size(), ""); - - const int ne = fesR->GetNE(); - const int nqe = ir->GetWeights().Size(); - - ElementTransformation *eltrans; - DofTransformation * doftrans; - const FiniteElement *fe = NULL; - Array vdofs; - - DenseMatrix trial_vshape; - - Vector Vx; - res.SetSize(rdim * rdim * rw.size()); - res = 0.0; - - // For the parallel case, we must get all DOFs of V on sampled elements. - CAROM::Matrix Vs; - - // Since V only has rows for true DOFs, we use a ParGridFunction to get all DOFs. - int eprev = -1; - int dof = 0; - int elemCount = 0; - - // First, find all sampled elements. - for (int i=0; iGetElementVDofs(e, vdofs); - if (dof > 0) - { - MFEM_VERIFY(dof == vdofs.Size(), "All elements must have same DOF size"); - } - dof = vdofs.Size(); - eprev = e; - elemCount++; - } - } - - // Now set Vs. - // Note that, ideally, these FOM data structures and operations should be - // avoided, but this only done in setup. - Vs.setSize(elemCount * dof, rdim); - ParGridFunction v_gf(fesR); - Vector vtrue(fesR->GetTrueVSize()); - - for (int j=0; jGetElementVDofs(e, vdofs); - - for (int k=0; k= 0) ? vdofs[k] : -1 - vdofs[k]; - const double vk = (vdofs[k] >= 0) ? v_gf[dofk] : -v_gf[dofk]; - Vs((elemCount * dof) + k, j) = vk; - } - - eprev = e; - elemCount++; - } - } - } - - eprev = -1; - elemCount = 0; - int spaceDim = 0; - - for (int i=0; iIntPoint(qpi); - - if (e != eprev) // Update element transformation - { - doftrans = fesR->GetElementVDofs(e, vdofs); - fe = fesR->GetFE(e); - eltrans = fesR->GetElementTransformation(e); - - if (doftrans) - { - MFEM_ABORT("TODO"); - } - - dof = fe->GetDof(); - spaceDim = eltrans->GetSpaceDim(); - trial_vshape.SetSize(dof, spaceDim); - Vx.SetSize(spaceDim); - - eprev = e; - elemCount++; - } - - // Integrate at the current point - - eltrans->SetIntPoint(&ip); - fe->CalcVShape(*eltrans, trial_vshape); - - double w = eltrans->Weight() * rw[i]; // using rw[i] instead of ip.weight - - // Note that the coefficient Q is omitted: w *= Q -> Eval(*eltrans, ip); - - for (int jx=0; jx const& rw, std::vector const& qp, const IntegrationRule *ir, Coefficient *Q, - CAROM::Matrix const& V, CAROM::Vector const& x, const int rank, Vector & res) -{ - const int rdim = V.numColumns(); - MFEM_VERIFY(rw.size() == qp.size(), ""); - MFEM_VERIFY(x.dim() == rdim, ""); - MFEM_VERIFY(V.numRows() == fesR->GetTrueVSize(), ""); - - MFEM_VERIFY(rank == 0, - "TODO: generalize to parallel. This uses full dofs in V, which has true dofs"); - - const int nqe = ir->GetWeights().Size(); - - ElementTransformation *eltrans; - DofTransformation * doftrans; - const FiniteElement *fe = NULL; - Array vdofs; - - DenseMatrix trial_vshape; - - Vector Vx; - res.SetSize(rdim); - res = 0.0; - - int eprev = -1; - int dof = 0; - int spaceDim = 0; - - // Note that the alternative version VectorFEMassIntegrator_ComputeReducedEQP_Fast - // of this function is optimized by storing some intermediate computations. - - for (int i=0; iIntPoint(qpi); - - if (e != eprev) // Update element transformation - { - doftrans = fesR->GetElementVDofs(e, vdofs); - fe = fesR->GetFE(e); - eltrans = fesR->GetElementTransformation(e); - - if (doftrans) - { - MFEM_ABORT("TODO"); - } - - dof = fe->GetDof(); - spaceDim = eltrans->GetSpaceDim(); - trial_vshape.SetSize(dof, spaceDim); - Vx.SetSize(spaceDim); - - eprev = e; - } - - // Integrate at the current point - - eltrans->SetIntPoint(&ip); - fe->CalcVShape(*eltrans, trial_vshape); - - double w = eltrans->Weight() * rw[i]; // using rw[i] instead of ip.weight - - if (Q) - { - w *= Q -> Eval(*eltrans, ip); - } - - // Lift Vx at ip. - Vx = 0.0; - for (int k=0; k= 0) ? vdofs[k] : -1 - vdofs[k]; - - double Vx_k = 0.0; - for (int j=0; j= 0) ? vdofs[l] : -1 - vdofs[l]; - const double s = (vdofs[l] >= 0) ? 1.0 : -1.0; - Vjk += s * V(dofl, j) * trial_vshape(l, k); - } - - rj += Vx[k] * Vjk; - } - - res[j] += w * rj; - } - } -} + CAROM::Matrix const& V, CAROM::Vector const& x, const int rank, Vector & res); void VectorFEMassIntegrator_ComputeReducedEQP_Fast(ParFiniteElementSpace *fesR, std::vector const& qp, const IntegrationRule *ir, Coefficient *Q, CAROM::Vector const& x, - Vector const& coef, Vector & res) -{ - const int rdim = x.dim(); - const int nqe = ir->GetWeights().Size(); - ElementTransformation *eltrans; - - res.SetSize(rdim); - res = 0.0; - - int eprev = -1; - - for (int i=0; iIntPoint(qpi); - - if (e != eprev) // Update element transformation - { - eltrans = fesR->GetElementTransformation(e); - eprev = e; - } - - eltrans->SetIntPoint(&ip); - const double q_ip = Q -> Eval(*eltrans, ip); - - for (int j=0; j const& rw, std::vector const& qp, const IntegrationRule *ir, CAROM::Matrix const& V, - Vector & res) -{ - // Assumption: fes is an L2 space, where true DOFs and full DOFs are the same. - // This allows for using V(dof,j), where dof is from fes->GetElementVDofs. - MFEM_VERIFY(fes->GetTrueVSize() == fes->GetTrueVSize(), ""); - - const int rdim = V.numColumns(); - MFEM_VERIFY(rw.size() == qp.size(), ""); - MFEM_VERIFY(V.numRows() == fes->GetTrueVSize(), ""); - - const int nqe = ir->GetWeights().Size(); - - ElementTransformation *eltrans; - DofTransformation * doftrans; - const FiniteElement *fe = NULL; - Array vdofs; - - Vector shape; - - int eprev = -1; - int dof = 0; - - for (int i=0; iIntPoint(qpi); - - if (e != eprev) // Update element transformation - { - doftrans = fes->GetElementVDofs(e, vdofs); - fe = fes->GetFE(e); - eltrans = fes->GetElementTransformation(e); - - if (doftrans) - { - MFEM_ABORT("TODO"); - } - - dof = fe->GetDof(); - - shape.SetSize(dof); - - eprev = e; - - if (i == 0) - { - res.SetSize(rdim * rw.size() * dof); - res = 0.0; - } - } - - // Integrate at the current point - - eltrans->SetIntPoint(&ip); - - fe->CalcPhysShape(*eltrans, shape); - - double w = eltrans->Weight() * rw[i]; // using rw[i] instead of ip.weight - - for (int j=0; j= 0) ? vdofs[l] : -1 - vdofs[l]; - const double s = (vdofs[l] >= 0) ? 1.0 : -1.0; - res[l + (j * dof) + (i * rdim * dof)] = w * s * V(dofl, - j) * shape(l) * shape(l); - } - } - } -} + Vector & res); // Based on MassIntegrator::AssembleElementMatrix void LinearMassIntegrator_ComputeReducedEQP(ParFiniteElementSpace *fes, std::vector const& rw, std::vector const& qp, const IntegrationRule *ir, Coefficient *Q, - CAROM::Matrix const& V, CAROM::Vector & res) -{ - const int rdim = V.numColumns(); - MFEM_VERIFY(rw.size() == qp.size(), ""); - MFEM_VERIFY(res.dim() == rdim, ""); - MFEM_VERIFY(V.numRows() == fes->GetTrueVSize(), ""); - - ParGridFunction f_gf(fes); - f_gf.ProjectCoefficient(*Q); - Vector sv(fes->GetTrueVSize()); - f_gf.GetTrueDofs(sv); - - const int nqe = ir->GetWeights().Size(); - - ElementTransformation *eltrans; - DofTransformation * doftrans; - const FiniteElement *fe = NULL; - Array vdofs; - - Vector shape; - - res = 0.0; - - int eprev = -1; - int dof = 0; - - // Note that the alternative version LinearMassIntegrator_ComputeReducedEQP_Fast - // of this function is optimized by storing some intermediate computations. - - for (int i=0; iIntPoint(qpi); - - if (e != eprev) // Update element transformation - { - doftrans = fes->GetElementVDofs(e, vdofs); - fe = fes->GetFE(e); - eltrans = fes->GetElementTransformation(e); - - if (doftrans) - { - MFEM_ABORT("TODO"); - } - - dof = fe->GetDof(); - - shape.SetSize(dof); - - eprev = e; - } - - // Integrate at the current point - - eltrans->SetIntPoint(&ip); - - fe->CalcPhysShape(*eltrans, shape); - - double w = eltrans->Weight() * rw[i]; // using rw[i] instead of ip.weight - - //w *= Q -> Eval(*eltrans, ip); - - for (int j=0; j= 0) ? vdofs[l] : -1 - vdofs[l]; - const double s = (vdofs[l] >= 0) ? 1.0 : -1.0; - //Vj += s * V(dofl, j) * shape(l); - Vj += s * V(dofl, j) * shape(l) * sv[dofl] * shape(l); - } - - res(j) += w * Vj; - } - } -} + CAROM::Matrix const& V, CAROM::Vector & res); void LinearMassIntegrator_ComputeReducedEQP_Fast(ParFiniteElementSpace *fes, ParGridFunction & f_gf, std::vector const& rw, std::vector const& qp, const IntegrationRule *ir, Coefficient *Q, CAROM::Matrix const& V, Vector const& coef, - CAROM::Vector & res) -{ - const int rdim = V.numColumns(); - MFEM_VERIFY(rw.size() == qp.size(), ""); - MFEM_VERIFY(V.numRows() == fes->GetTrueVSize(), ""); - - const int nqe = ir->GetWeights().Size(); - - DofTransformation * doftrans; - const FiniteElement *fe = NULL; - Array vdofs; - - res = 0.0; - - int eprev = -1; - int dof = 0; - - for (int i=0; iIntPoint(qpi); - - if (e != eprev) // Update element transformation - { - doftrans = fes->GetElementVDofs(e, vdofs); - fe = fes->GetFE(e); - - if (doftrans) - { - MFEM_ABORT("TODO"); - } - - dof = fe->GetDof(); - - eprev = e; - - f_gf.ProjectCoefficient(*Q, vdofs); - - if (i == 0) - { - MFEM_VERIFY(coef.Size() == rdim * rw.size() * dof, ""); - } - } - - // Integrate at the current point - - for (int j=0; j= 0) ? vdofs[l] : -1 - vdofs[l]; - Vj += f_gf[dofl] * coef[l + (j * dof) + (i * rdim * dof)]; - } - - res(j) += Vj; - } - } -} + CAROM::Vector & res); // Compute a(p) u . v at all quadrature points on the given element. Coefficient Q is a(p). void ComputeElementRowOfG(const IntegrationRule *ir, Array const& vdofs, Coefficient *Q, Vector const& u, Vector const& v, - FiniteElement const& fe, ElementTransformation & Trans, Vector & r) -{ - MFEM_VERIFY(r.Size() == ir->GetNPoints(), ""); - int dof = fe.GetDof(); - int spaceDim = Trans.GetSpaceDim(); - - Vector u_i(spaceDim); - Vector v_i(spaceDim); - - DenseMatrix trial_vshape(dof, spaceDim); - - for (int i = 0; i < ir->GetNPoints(); i++) - { - const IntegrationPoint &ip = ir->IntPoint(i); - - Trans.SetIntPoint(&ip); - - fe.CalcVShape(Trans, trial_vshape); - - double w = Trans.Weight(); - - u_i = 0.0; - v_i = 0.0; - - for (int j=0; j= 0) ? vdofs[j] : -1 - vdofs[j]; - const double s = (vdofs[j] >= 0) ? 1.0 : -1.0; - for (int k=0; k Eval (Trans, ip); - } - - r[i] = 0.0; - for (int k=0; k const& vdofs, Vector const& s, Vector const& p, - FiniteElement const& fe, ElementTransformation & Trans, Vector & r) -{ - MFEM_VERIFY(r.Size() == ir->GetNPoints(), ""); - int dof = fe.GetDof(); - - Vector shape(dof); - - for (int i = 0; i < ir->GetNPoints(); i++) - { - const IntegrationPoint &ip = ir->IntPoint(i); - - Trans.SetIntPoint(&ip); - - fe.CalcPhysShape(Trans, shape); - - double w = Trans.Weight(); - - double u_i = 0.0; - double v_i = 0.0; - - for (int j=0; j= 0) ? vdofs[j] : -1 - vdofs[j]; - const double sd = (vdofs[j] >= 0) ? 1.0 : -1.0; - u_i += sd * p[dofj] * shape(j); - v_i += sd * s[dofj] * shape(j); - } - - r[i] = w * u_i * v_i; - } -} + FiniteElement const& fe, ElementTransformation & Trans, Vector & r); // Precondition Gt by (V^T M V)^{-1} (of size NB x NB). void PreconditionNNLS(ParFiniteElementSpace *fespace, BilinearFormIntegrator *massInteg, const CAROM::Matrix* B, const int snapshot, - CAROM::Matrix & Gt) -{ - const int NB = B->numColumns(); - const int NQ = Gt.numRows(); - const int nsnap = Gt.numColumns() / NB; - - ParBilinearForm M(fespace); - - M.AddDomainIntegrator(massInteg); - M.Assemble(); - M.Finalize(); - HypreParMatrix *Mmat = M.ParallelAssemble(); - - CAROM::Matrix Mhat(NB, NB, false); - ComputeCtAB(*Mmat, *B, *B, Mhat); - Mhat.inverse(); - - CAROM::Vector PG(NB, false); - - for (int i = (snapshot >= 0 ? snapshot : 0); - i < (snapshot >= 0 ? snapshot+1 : nsnap); ++i) - { - for (int m=0; mGetNPoints(); - const int ne = fespace_R->GetNE(); - const int NB = BR->numColumns(); - const int NQ = ne * nqe; - const int nsnap = BR_snapshots->numColumns(); - - MFEM_VERIFY(nsnap == BW_snapshots->numColumns() || - nsnap + nsets == BW_snapshots->numColumns(), ""); - MFEM_VERIFY(BR->numRows() == BR_snapshots->numRows(), ""); - MFEM_VERIFY(BR->numRows() == fespace_R->GetTrueVSize(), ""); - MFEM_VERIFY(BW_snapshots->numRows() == fespace_W->GetTrueVSize(), ""); - - const bool skipFirstW = (nsnap + nsets == BW_snapshots->numColumns()); - - // Compute G of size (NB * nsnap) x NQ, but only store its transpose Gt. - CAROM::Matrix Gt(NQ, NB * nsnap, true); - - // For 0 <= j < NB, 0 <= i < nsnap, 0 <= e < ne, 0 <= m < nqe, - // G(j + (i*NB), (e*nqe) + m) - // is the coefficient of v_j^T M(p_i) V 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 p_i(BW_snapshots->numRows()); - Vector v_i(BR_snapshots->numRows()); - Vector v_j(BR->numRows()); - - Vector r(nqe); - - int skip = 0; - const int nsnapPerSet = nsnap / nsets; - if (skipFirstW) - { - MFEM_VERIFY(nsets * nsnapPerSet == nsnap, ""); - skip = 1; - } - - for (int i=0; inumRows(); ++j) - p_i[j] = (*BW_snapshots)(j, i + skip); - - for (int j = 0; j < BR_snapshots->numRows(); ++j) - v_i[j] = (*BR_snapshots)(j, i); - - if (skipFirstW && i > 0 && i % nsnapPerSet == 0) - skip++; - - // Set grid function for a(p) - ParGridFunction p_gf(fespace_W); - - p_gf.SetFromTrueDofs(p_i); - - GridFunctionCoefficient p_coeff(&p_gf); - TransformedCoefficient a_coeff(&p_coeff, NonlinearCoefficient); - - ParGridFunction vi_gf(fespace_R); - vi_gf.SetFromTrueDofs(v_i); - - for (int j=0; jnumRows(); ++k) - v_j[k] = (*BR)(k, j); - - ParGridFunction vj_gf(fespace_R); - vj_gf.SetFromTrueDofs(v_j); - - // TODO: is it better to make the element loop the outer loop? - for (int e=0; e vdofs; - DofTransformation *doftrans = fespace_R->GetElementVDofs(e, vdofs); - const FiniteElement &fe = *fespace_R->GetFE(e); - ElementTransformation *eltrans = fespace_R->GetElementTransformation(e); - - ComputeElementRowOfG(ir0, vdofs, &a_coeff, vi_gf, vj_gf, fe, *eltrans, r); - - for (int m=0; m const& w_el = ir0->GetWeights(); - MFEM_VERIFY(w_el.Size() == nqe, ""); - - CAROM::Vector w(ne * nqe, true); - - for (int i=0; iGetNPoints(); - const int ne = fespace_W->GetNE(); - const int NB = BW->numColumns(); - const int NQ = ne * nqe; - const int nsnap = BS_snapshots->numColumns(); - const int nrows = BW->numRows(); - - MFEM_VERIFY(BS_snapshots->numRows() == BW->numRows(), ""); - - // Compute G of size (NB * nsnap) x NQ, but only store its transpose Gt. - CAROM::Matrix Gt(NQ, NB * nsnap, true); - - Vector s_i(nrows); - Vector p_j(nrows); - - Vector r(nqe); - - for (int i=0; inumRows(); ++j) - s_i[j] = (*BS_snapshots)(j, i); - - for (int j=0; jnumRows(); ++k) - p_j[k] = (*BW)(k, j); - - // TODO: is it better to make the element loop the outer loop? - for (int e=0; e vdofs; - DofTransformation *doftrans = fespace_W->GetElementVDofs(e, vdofs); - const FiniteElement &fe = *fespace_W->GetFE(e); - ElementTransformation *eltrans = fespace_W->GetElementTransformation(e); - - ComputeElementRowOfG_Source(ir0, vdofs, s_i, p_j, fe, *eltrans, r); - - for (int m=0; m const& w_el = ir0->GetWeights(); - MFEM_VERIFY(w_el.Size() == nqe, ""); - - CAROM::Vector w(ne * nqe, true); - - for (int i=0; i elements; - Vector elemCount(pmesh->GetNE()); - elemCount = 0.0; - - for (int i=0; i 1.0e-12) - { - const int e = i / nqe; // Element index - elements.insert(e); - elemCount[e] += 1.0; - } - } + CAROM::Vector const& eqpSol); - // Empty sets, since EQP on samples inside elements. - std::set faces, edges, vertices; - CAROM::SampleVisualization(pmesh, elements, elements, faces, edges, - vertices, "EQPvis", &elemCount); -} +#endif diff --git a/examples/prom/nonlinear_elasticity_global_rom_eqp.cpp b/examples/prom/nonlinear_elasticity_global_rom_eqp.cpp new file mode 100644 index 000000000..077ac3377 --- /dev/null +++ b/examples/prom/nonlinear_elasticity_global_rom_eqp.cpp @@ -0,0 +1,811 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// Functions used by mixed_nonlinear_diffusion.cpp with EQP. +#include "mfem/Utilities.hpp" +// #include "fem/nonlininteg.hpp" +#include +using namespace mfem; +using namespace std; + +#include "nonlinear_elasticity_global_rom_eqp.hpp" + +class HyperelasticOperator : public TimeDependentOperator +{ + +protected: + ParBilinearForm *M, *S; + + CGSolver M_solver; // Krylov solver for inverting the mass matrix M + HypreSmoother M_prec; // Preconditioner for the mass matrix M + +public: + HyperelasticOperator(ParFiniteElementSpace &f, Array &ess_tdof_list_, + double visc, double mu, double K); + + /// Compute the right-hand side of the ODE system. + virtual void Mult(const Vector &vx, Vector &dvx_dt) const; + + double ElasticEnergy(const ParGridFunction &x) const; + double KineticEnergy(const ParGridFunction &v) const; + void GetElasticEnergyDensity(const ParGridFunction &x, + ParGridFunction &w) const; + + mutable Vector H_sp; + mutable Vector dvxdt_sp; + + ParFiniteElementSpace &fespace; + double viscosity; + Array ess_tdof_list; + ParNonlinearForm *H; + HyperelasticModel *model; + mutable Vector z; // auxiliary vector + mutable Vector z2; // auxiliary vector + HypreParMatrix *Mmat; // Mass matrix from ParallelAssemble() + HypreParMatrix Smat; + + virtual ~HyperelasticOperator(); +}; + +void GetEQPCoefficients_HyperelasticNLFIntegrator(ParFiniteElementSpace *fesR, + std::vector const &rw, std::vector const &qp, + const IntegrationRule *ir, NeoHookeanModel *model, + CAROM::Matrix const &V_v, const int rank, Vector &coef, Vector &DS_coef) +{ + const int rvdim = V_v.numColumns(); + const int fomdim = V_v.numRows(); + MFEM_VERIFY(rw.size() == qp.size(), ""); + MFEM_VERIFY(fomdim == fesR->GetTrueVSize(), ""); + + MFEM_VERIFY(rank == 0, + "TODO: generalize to parallel. This uses full dofs in V, which has true dofs"); + + const int nqe = ir->GetWeights().Size(); + + DofTransformation *doftrans; + ElementTransformation *eltrans; + Array vdofs; + + int eprev = -1; + int dof = 0; + int dim = 0; + + // Get prolongation matrix + const Operator *P = fesR->GetProlongationMatrix(); + + // Vector to be prolongated + Vector vj(fomdim); + + // Prolongated vector + Vector p_vj(P->Height()); + + // Element vector + Vector vj_e; + // Get the element vector size + doftrans = fesR->GetElementVDofs(0, vdofs); + int elvect_size = vdofs.Size(); + + // Coefficient vector + coef.SetSize(elvect_size * rw.size() * rvdim); + coef = 0.0; + + // Vector for storing DS + const FiniteElement *fe = fesR->GetFE(0); + dof = fe->GetDof(); + dim = fe->GetDim(); + DenseMatrix DSh(dof, dim); + DenseMatrix DS(dof, dim); + DenseMatrix Jrt(dim); + DS_coef.SetSize(dof * dim * rw.size() * rvdim); + DS_coef = 0.0; + int index = 0; + + // For every basis vector + for (int j = 0; j < rvdim; ++j) + { + // Get basis vector and prolongate + for (int k = 0; k < V_v.numRows(); ++k) + vj[k] = V_v(k, j); + P->Mult(vj, p_vj); + + eprev = -1; + + // For every quadrature weight + for (int i = 0; i < rw.size(); ++i) // NOTE: i < 9 + { + const int e = qp[i] / nqe; // Element index + // Local (element) index of the quadrature point + const int qpi = qp[i] - (e * nqe); + const IntegrationPoint &ip = ir->IntPoint(qpi); + + if (e != eprev) // Update element transformation + { + doftrans = fesR->GetElementVDofs(e, vdofs); + eltrans = fesR->GetElementTransformation(e); + + if (doftrans) + { + MFEM_ABORT("TODO"); + } + + // Get element vectors + p_vj.GetSubVector(vdofs, vj_e); + eprev = e; + } + + // Set integration point in the element transformation + eltrans->SetIntPoint(&ip); + model->SetTransformation(*eltrans); + // Get the transformation weight + double t = eltrans->Weight(); + // Calculate r[i] = ve_j^T * elvect + for (int k = 0; k < elvect_size; k++) + { + coef[k + (i * elvect_size) + (j * rw.size() * elvect_size)] = vj_e[k] * rw[i] * t; + } + + // Calculate DS and store + CalcInverse(eltrans->Jacobian(), Jrt); + fe->CalcDShape(ip, DSh); + Mult(DSh, Jrt, DS); + for (int ii = 0; ii < dof; ++ii) + { + for (int jj = 0; jj < dim; ++jj) + { + index = jj + ii * dim; + DS_coef[index + (i * dof*dim) + (j * rw.size() * dof*dim)] = DS.Elem(ii,jj); + } + } + } + } +} + +void HyperelasticNLFIntegrator_ComputeReducedEQP(ParFiniteElementSpace *fesR, + std::vector const &rw, std::vector const &qp, + const IntegrationRule *ir, NeoHookeanModel *model, const Vector *x0, + CAROM::Matrix const &V_x, CAROM::Matrix const &V_v, CAROM::Vector const &x, const int rank, Vector &res) +{ + const int rxdim = V_x.numColumns(); + const int rvdim = V_v.numColumns(); + const int fomdim = V_v.numRows(); + MFEM_VERIFY(rw.size() == qp.size(), ""); + MFEM_VERIFY(x.dim() == rxdim, ""); + MFEM_VERIFY(V_x.numRows() == fesR->GetTrueVSize(), ""); + + MFEM_VERIFY(rank == 0, + "TODO: generalize to parallel. This uses full dofs in V, which has true dofs"); + + const int nqe = ir->GetWeights().Size(); + + ElementTransformation *eltrans; + DofTransformation *doftrans; + const FiniteElement *fe = NULL; + Array vdofs; + + res.SetSize(rxdim); + res = 0.0; + + int eprev = -1; + int dof = 0; + int dim = 0; + + // Get prolongation matrix + const Operator *P = fesR->GetProlongationMatrix(); + + // Vectors to be prolongated + CAROM::Vector *Vx_librom_temp = new CAROM::Vector(fomdim, true); + Vector *Vx_temp = new Vector(&((*Vx_librom_temp)(0)), fomdim); + Vector Vx(fomdim); + Vector vj(fomdim); + + // Prolongated vectors + Vector p_Vx(P->Height()); + Vector p_vj(P->Height()); + + // Element vectors + Vector Vx_e; + Vector vj_e; + + // Lift x, add x0 and prolongate result + V_x.mult(x, Vx_librom_temp); + add(*Vx_temp, *x0, Vx); + P->Mult(Vx, p_Vx); + + // Initialize nonlinear operator storage + // Assuming all elements have the same dim and n_dof + fe = fesR->GetFE(0); + dof = fe->GetDof(); + dim = fe->GetDim(); + DenseMatrix DSh(dof, dim); + DenseMatrix DS(dof, dim); + DenseMatrix Jrt(dim); + DenseMatrix Jpt(dim); + DenseMatrix P_f(dim); + DenseMatrix PMatI; // Extract element dofs + DenseMatrix PMatO; + Vector elvect(dof * dim); + PMatO.UseExternalData(elvect.GetData(), dof, dim); + + // For every basis vector + for (int j = 0; j < rvdim; ++j) + + { + // Get basis vector and prolongate + for (int k = 0; k < V_v.numRows(); ++k) + vj[k] = V_v(k, j); + P->Mult(vj, p_vj); + res[j] = 0.0; + + eprev = -1; + + // For every quadrature weight + for (int i = 0; i < rw.size(); ++i) // NOTE: i < 9 + { + const int e = qp[i] / nqe; // Element index + // Local (element) index of the quadrature point + const int qpi = qp[i] - (e * nqe); + const IntegrationPoint &ip = ir->IntPoint(qpi); + + if (e != eprev) // Update element transformation + { + doftrans = fesR->GetElementVDofs(e, vdofs); + fe = fesR->GetFE(e); + eltrans = fesR->GetElementTransformation(e); + + dof = fe->GetDof(); // Get number of dofs in element + dim = fe->GetDim(); + + if (doftrans) + { + MFEM_ABORT("TODO"); + } + + // Get element vectors + p_Vx.GetSubVector(vdofs, Vx_e); + p_vj.GetSubVector(vdofs, vj_e); + eprev = e; + } + + // Integration at ip + PMatI.UseExternalData(Vx_e.GetData(), dof, dim); + + elvect = 0.0; + + // Set integration point in the element transformation + eltrans->SetIntPoint(&ip); + model->SetTransformation(*eltrans); + + // Get the transformation weight + double t = eltrans->Weight(); + + // Compute action of nonlinear operator + CalcInverse(eltrans->Jacobian(), Jrt); + fe->CalcDShape(ip, DSh); + Mult(DSh, Jrt, DS); + MultAtB(PMatI, DS, Jpt); + model->EvalP(Jpt, P_f); + P_f *= (t * rw[i]); // NB: Not by ip.weight + AddMultABt(DS, P_f, PMatO); + + // Calculate r[i] = ve_j^T * elvect + for (int k = 0; k < elvect.Size(); k++) + { + res[j] += vj_e[k] * elvect[k]; + } + } + } +} + +void HyperelasticNLFIntegrator_ComputeReducedEQP_Fast(ParFiniteElementSpace *fesR, + std::vector const &rw, std::vector const &qp, const IntegrationRule *ir, NeoHookeanModel *model, + const Vector *x0, CAROM::Matrix const &V_x, CAROM::Matrix const &V_v, CAROM::Vector const &x, + Vector const &coef, Vector const &DS_coef, const int rank, Vector &res) +{ + + const int rxdim = V_x.numColumns(); + const int rvdim = V_v.numColumns(); + const int fomdim = V_x.numRows(); + MFEM_VERIFY(x.dim() == rxdim, ""); + MFEM_VERIFY(V_x.numRows() == fesR->GetTrueVSize(), ""); + + MFEM_VERIFY(rank == 0, + "TODO: generalize to parallel. This uses full dofs in V, which has true dofs"); + + const int nqe = ir->GetWeights().Size(); + + ElementTransformation *eltrans; + DofTransformation *doftrans; + const FiniteElement *fe = NULL; + Array vdofs; + + res.SetSize(rxdim); + res = 0.0; + + int eprev = -1; + int dof = 0; + int dim = 0; + + // Get prolongation matrix + const Operator *P = fesR->GetProlongationMatrix(); + + // Vectors to be prolongated + CAROM::Vector *Vx_librom_temp = new CAROM::Vector(fomdim, true); + Vector *Vx_temp = new Vector(&((*Vx_librom_temp)(0)), fomdim); + Vector Vx(fomdim); + + // Prolongated vectors + Vector p_Vx(P->Height()); + + // Element vectors + Vector Vx_e; + + // Lift x, add x0 and prolongate result + V_x.mult(x, Vx_librom_temp); + add(*Vx_temp, *x0, Vx); + P->Mult(Vx, p_Vx); + + // Initialize nonlinear operator storage + // Assuming all elements have the same dim and n_dof + fe = fesR->GetFE(0); + dof = fe->GetDof(); + dim = fe->GetDim(); + int index = 0; + DenseMatrix DS(dof, dim); + DenseMatrix Jpt(dim); + DenseMatrix P_f(dim); + DenseMatrix PMatI; // Extract element dofs + DenseMatrix PMatO; + Vector elvect(dof * dim); + PMatO.UseExternalData(elvect.GetData(), dof, dim); + + // For every basis vector + for (int j = 0; j < rvdim; ++j) + + { + + eprev = -1; + + // For every quadrature weight + for (int i = 0; i < qp.size(); ++i) // NOTE: i < 9 + { + const int e = qp[i] / nqe; // Element index + // Local (element) index of the quadrature point + const int qpi = qp[i] - (e * nqe); + const IntegrationPoint &ip = ir->IntPoint(qpi); + + if (e != eprev) // Update element transformation + { + doftrans = fesR->GetElementVDofs(e, vdofs); + fe = fesR->GetFE(e); + eltrans = fesR->GetElementTransformation(e); + + dof = fe->GetDof(); // Get number of dofs in element + dim = fe->GetDim(); + + if (doftrans) + { + MFEM_ABORT("TODO"); + } + + // Get element vectors + p_Vx.GetSubVector(vdofs, Vx_e); + eprev = e; + } + + // Integration at ip + elvect = 0.0; + PMatI.UseExternalData(Vx_e.GetData(), dof, dim); + + // Set integration point in the element transformation + eltrans->SetIntPoint(&ip); + model->SetTransformation(*eltrans); + + for (int ii = 0; ii < dof; ++ii) + { + for (int jj = 0; jj < dim; ++jj) + { + index = jj + ii * dim; + DS.Elem(ii,jj) = DS_coef[index + (i * elvect.Size()) + (j * qp.size() * elvect.Size())]; + } + } + + MultAtB(PMatI, DS, Jpt); + model->EvalP(Jpt, P_f); + AddMultABt(DS, P_f, PMatO); + + // Calculate r[i] = ve_j^T * elvect + // coef is size len(vdofs) * rvdim * rw.size + for (int k = 0; k < elvect.Size(); k++) + { + res[j] += coef[k + (i * elvect.Size()) + (j * qp.size() * elvect.Size())] * elvect[k]; + } + } + } +} + +bool fileExists(const std::string &filename) +{ + std::ifstream file(filename); + return file.good(); +} + +void save_CAROM_Vector(const CAROM::Vector &vector, const std::string &filename) +{ + std::ofstream file(filename); + if (file.is_open()) + { + for (int i = 0; i < vector.dim(); ++i) + { + file << vector(i) << "\n"; + } + file.close(); + std::cout << "Vector saved as " << filename << " successfully." << std::endl; + } + else + { + std::cerr << "Error: Unable to open file " << filename << " for writing." << std::endl; + } +} + +void ComputeElementRowOfG(const IntegrationRule *ir, Array const &vdofs, + Vector const &ve_j, NeoHookeanModel *model, Vector const &elfun, + FiniteElement const &fe, ElementTransformation &Trans, Vector &r) +{ + MFEM_VERIFY(r.Size() == ir->GetNPoints(), ""); + int dof = fe.GetDof(); // Get number of dofs in element + int dim = fe.GetDim(); + + // Initialize nonlinear operator matrices (there is probably a better way) + DenseMatrix DSh(dof, dim); + DenseMatrix DS(dof, dim); + DenseMatrix Jrt(dim); + DenseMatrix Jpt(dim); + DenseMatrix P(dim); + DenseMatrix PMatI; // Extract element dofs + PMatI.UseExternalData(elfun.GetData(), dof, dim); + DenseMatrix PMatO; + Vector elvect(dof * dim); + PMatO.UseExternalData(elvect.GetData(), dof, dim); + + model->SetTransformation(Trans); + + // For each integration point + for (int i = 0; i < ir->GetNPoints(); i++) + { + elvect = 0.0; + // Get integration point + const IntegrationPoint &ip = ir->IntPoint(i); + + // Set integration point in the element transformation + Trans.SetIntPoint(&ip); + + // Get the transformation weight + double t = Trans.Weight(); + + // Compute action of nonlinear operator + CalcInverse(Trans.Jacobian(), Jrt); + fe.CalcDShape(ip, DSh); + Mult(DSh, Jrt, DS); + MultAtB(PMatI, DS, Jpt); + model->EvalP(Jpt, P); + P *= t; // NB: Not by ip.weight + AddMultABt(DS, P, PMatO); + + r[i] = 0.0; + + // Calculate r[i] = ve_j^T * elvect + for (int k = 0; k < elvect.Size(); k++) + { + r[i] += ve_j[k] * elvect[k]; + } + } +} + +void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, + CAROM::Vector const &w, CAROM::Matrix &Gt, + CAROM::Vector &sol) +{ + 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); + + CAROM::Vector rhs_lb(rhs_ub); + CAROM::Vector rhs_Gw(rhs_ub); + + const double delta = 1.0e-11; + for (int i = 0; i < rhs_ub.dim(); ++i) + { + rhs_lb(i) -= delta; + rhs_ub(i) += delta; + } + + // nnls.normalize_constraints(Gt, rhs_lb, rhs_ub); + nnls.solve_parallel_with_scalapack(Gt, rhs_lb, rhs_ub, sol); + + int nnz = 0; + for (int i = 0; i < sol.dim(); ++i) + { + if (sol(i) != 0.0) + { + nnz++; + } + } + + cout << rank << ": Number of nonzeros in NNLS solution: " << nnz + << ", out of " << sol.dim() << endl; + + MPI_Allreduce(MPI_IN_PLACE, &nnz, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + + if (rank == 0) + cout << "Global number of nonzeros in NNLS solution: " << nnz << endl; + + // Check residual of NNLS solution + CAROM::Vector res(Gt.numColumns(), false); + Gt.transposeMult(sol, res); + + const double normGsol = res.norm(); + const double normRHS = rhs_Gw.norm(); + + res -= rhs_Gw; + const double relNorm = res.norm() / std::max(normGsol, normRHS); + cout << rank << ": relative residual norm for NNLS solution of Gs = Gw: " << relNorm << endl; +} + +void SetupEQP_snapshots(const IntegrationRule *ir0, const int rank, + ParFiniteElementSpace *fespace_X, + const int nsets, const CAROM::Matrix *BV, + const CAROM::Matrix *BX_snapshots, + const bool precondition, const double nnls_tol, + const int maxNNLSnnz, NeoHookeanModel *model, + CAROM::Vector &sol, CAROM::Vector *window_ids) +{ + const int nqe = ir0->GetNPoints(); + const int ne = fespace_X->GetNE(); + const int NB = BV->numColumns(); + const int NQ = ne * nqe; + const int nsnap = BX_snapshots->numColumns(); + + MFEM_VERIFY(nsnap == BX_snapshots->numColumns() || + nsnap + nsets == BX_snapshots->numColumns(), // Q: nsets? + ""); + MFEM_VERIFY(BV->numRows() == BX_snapshots->numRows(), ""); + MFEM_VERIFY(BV->numRows() == fespace_X->GetTrueVSize(), ""); + + const bool skipFirstW = (nsnap + nsets == BX_snapshots->numColumns()); + + // Compute G of size (NB * nsnap) x NQ, but only store its transpose Gt. + CAROM::Matrix Gt(NQ, NB * nsnap, true); + int current_size = 0; + int i_start = 0; + int outer_loop_length = 0; + if (!window_ids) + { + current_size = nsnap; + outer_loop_length = 1; + } + else + { + outer_loop_length = window_ids->dim() - 1; + } + + // For 0 <= j < NB, 0 <= i < nsnap, 0 <= e < ne, 0 <= m < nqe, + // G(j + (i*NB), (e*nqe) + m) + // is the coefficient of v_j^T H(x_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 x_i(BX_snapshots->numRows()); + Vector v_j(BV->numRows()); + + Vector r(nqe); + + int skip = 0; + const int nsnapPerSet = nsnap / nsets; + if (skipFirstW) + { + MFEM_VERIFY(nsets * nsnapPerSet == nsnap, ""); + skip = 1; + } + + // Get prolongation matrix + const Operator *P = fespace_X->GetProlongationMatrix(); + if (!P) + { + MFEM_ABORT("P is null, generalize to serial case") + } + + // Outer loop for time windowing + for (int oi = 0; oi < outer_loop_length; ++oi) + { + if (window_ids) + { + i_start = window_ids->item(oi); + current_size = window_ids->item(oi + 1) - window_ids->item(oi) + 1; + cout << "i_start = " << i_start << endl; + Gt.setSize(NQ, NB * current_size); + } + + // For every snapshot in batch + for (int i = i_start; i < current_size; ++i) + { + // Set the sampled dofs from the snapshot matrix + for (int j = 0; j < BX_snapshots->numRows(); ++j) + x_i[j] = (*BX_snapshots)(j, i); + + // Get prolongated dofs + Vector px_i; + Vector elfun; + + px_i.SetSize(P->Height()); + P->Mult(x_i, px_i); + + if (skipFirstW && i > 0 && i % nsnapPerSet == 0) + skip++; + + // For each basis vector + for (int j = 0; j < NB; ++j) + { + + // Get basis vector + for (int k = 0; k < BV->numRows(); ++k) + v_j[k] = (*BV)(k, j); + + // Get prolongated dofs + Vector pv_j; + Vector ve_j; + + pv_j.SetSize(P->Height()); + P->Mult(v_j, pv_j); + + // TODO: is it better to make the element loop the outer loop? + // For each element + for (int e = 0; e < ne; ++e) + { + // Get element and its dofs and transformation. + Array vdofs; + DofTransformation *doftrans = fespace_X->GetElementVDofs(e, vdofs); + const FiniteElement &fe = *fespace_X->GetFE(e); + ElementTransformation *eltrans = fespace_X->GetElementTransformation(e); + px_i.GetSubVector(vdofs, elfun); + pv_j.GetSubVector(vdofs, ve_j); + if (doftrans) + { + MFEM_ABORT("Doftrans is true, make corresponding edits") + } + + // Compute the row of G corresponding to element e, store in r + ComputeElementRowOfG(ir0, vdofs, ve_j, model, elfun, fe, *eltrans, r); + + for (int m = 0; m < nqe; ++m) + Gt((e * nqe) + m, j + (i * NB)) = r[m]; + } + } + + if (precondition) + { + MFEM_ABORT("TODO"); + } + } // Loop (i) over snapshots + + Array const &w_el = ir0->GetWeights(); + MFEM_VERIFY(w_el.Size() == nqe, ""); + + CAROM::Vector w(ne * nqe, true); + + for (int i = 0; i < ne; ++i) + { + for (int j = 0; j < nqe; ++j) + w((i * nqe) + j) = w_el[j]; + } + + SolveNNLS(rank, nnls_tol, maxNNLSnnz, w, Gt, sol); + if (window_ids) + { + save_CAROM_Vector(sol, "sol_" + std::to_string(oi) + ".csv"); + } + } +} + +void WriteMeshEQP(ParMesh *pmesh, const int myid, const int nqe, + CAROM::Vector const &eqpSol) +{ + // Find the elements with quadrature points in eqpSol. + std::set elements; + Vector elemCount(pmesh->GetNE()); + elemCount = 0.0; + + for (int i = 0; i < eqpSol.dim(); ++i) + { + if (eqpSol(i) > 1.0e-12) + { + const int e = i / nqe; // Element index + elements.insert(e); + elemCount[e] += 1.0; + } + } + + // Empty sets, since EQP on samples inside elements. + std::set faces, edges, vertices; + CAROM::SampleVisualization(pmesh, elements, elements, faces, edges, + vertices, "EQPvis", &elemCount); +} + +void get_window_ids(int n_step, int n_window, CAROM::Vector *ids) +{ + int window_size = std::round(n_step / n_window); + + int res = n_step - (n_window * (window_size - 1) + 1); + int res_up = res - n_window; + int ctr = -1; + + for (int i = 1; i <= n_window + 1; ++i) + { + ids->item(i - 1) = (i - 1) * window_size + 1 - (i - 1); + if (res_up > 0) + { + if (i <= res - res_up) + { + ctr += 1; + } + if (i > 2 && i <= res_up + 1) + { + ctr += 1; + } + } + else + { + if (i <= res + 1) + { + ctr += 1; + } + } + ids->item(i - 1) += ctr; + } +} + +void load_CAROM_vector(const std::string &filename, CAROM::Vector &vector) +{ + std::ifstream file(filename); + std::string line; + std::vector data; + + while (std::getline(file, line)) + { + std::stringstream ss(line); + std::string value; + + while (std::getline(ss, value, ',')) + { + data.push_back(std::stod(value)); + } + } + + // Set the size of the vector + int size = data.size(); + vector.setSize(size); + + // Copy the data into the vector + for (int i = 0; i < size; ++i) + { + vector(i) = data[i]; + } +} + +void get_EQPsol(const int current_window, CAROM::Vector *load_eqpsol) +{ + string filename = "sol_" + std::to_string(current_window) + ".csv"; + if (fileExists(filename)) + { + load_CAROM_vector(filename, *load_eqpsol); + } +} diff --git a/examples/prom/nonlinear_elasticity_global_rom_eqp.hpp b/examples/prom/nonlinear_elasticity_global_rom_eqp.hpp index bc6727485..d50db71c0 100644 --- a/examples/prom/nonlinear_elasticity_global_rom_eqp.hpp +++ b/examples/prom/nonlinear_elasticity_global_rom_eqp.hpp @@ -8,6 +8,8 @@ * *****************************************************************************/ +#ifndef NLEL_EQP +#define NLEL_EQP // Functions used by mixed_nonlinear_diffusion.cpp with EQP. #include "mfem/Utilities.hpp" // #include "fem/nonlininteg.hpp" @@ -15,6 +17,7 @@ using namespace mfem; using namespace std; +//TODO: Check what to do with class definition class HyperelasticOperator : public TimeDependentOperator { @@ -57,509 +60,29 @@ class HyperelasticOperator : public TimeDependentOperator void GetEQPCoefficients_HyperelasticNLFIntegrator(ParFiniteElementSpace *fesR, std::vector const &rw, std::vector const &qp, const IntegrationRule *ir, NeoHookeanModel *model, - CAROM::Matrix const &V_v, const int rank, Vector &coef, Vector &DS_coef) -{ - const int rvdim = V_v.numColumns(); - const int fomdim = V_v.numRows(); - MFEM_VERIFY(rw.size() == qp.size(), ""); - MFEM_VERIFY(fomdim == fesR->GetTrueVSize(), ""); - - MFEM_VERIFY(rank == 0, - "TODO: generalize to parallel. This uses full dofs in V, which has true dofs"); - - const int nqe = ir->GetWeights().Size(); - - DofTransformation *doftrans; - ElementTransformation *eltrans; - Array vdofs; - - int eprev = -1; - int dof = 0; - int dim = 0; - - // Get prolongation matrix - const Operator *P = fesR->GetProlongationMatrix(); - - // Vector to be prolongated - Vector vj(fomdim); - - // Prolongated vector - Vector p_vj(P->Height()); - - // Element vector - Vector vj_e; - // Get the element vector size - doftrans = fesR->GetElementVDofs(0, vdofs); - int elvect_size = vdofs.Size(); - - // Coefficient vector - coef.SetSize(elvect_size * rw.size() * rvdim); - coef = 0.0; - - // Vector for storing DS - const FiniteElement *fe = fesR->GetFE(0); - dof = fe->GetDof(); - dim = fe->GetDim(); - DenseMatrix DSh(dof, dim); - DenseMatrix DS(dof, dim); - DenseMatrix Jrt(dim); - DS_coef.SetSize(dof * dim * rw.size() * rvdim); - DS_coef = 0.0; - int index = 0; - - // For every basis vector - for (int j = 0; j < rvdim; ++j) - { - // Get basis vector and prolongate - for (int k = 0; k < V_v.numRows(); ++k) - vj[k] = V_v(k, j); - P->Mult(vj, p_vj); - - eprev = -1; - - // For every quadrature weight - for (int i = 0; i < rw.size(); ++i) // NOTE: i < 9 - { - const int e = qp[i] / nqe; // Element index - // Local (element) index of the quadrature point - const int qpi = qp[i] - (e * nqe); - const IntegrationPoint &ip = ir->IntPoint(qpi); - - if (e != eprev) // Update element transformation - { - doftrans = fesR->GetElementVDofs(e, vdofs); - eltrans = fesR->GetElementTransformation(e); - - if (doftrans) - { - MFEM_ABORT("TODO"); - } - - // Get element vectors - p_vj.GetSubVector(vdofs, vj_e); - eprev = e; - } - - // Set integration point in the element transformation - eltrans->SetIntPoint(&ip); - model->SetTransformation(*eltrans); - // Get the transformation weight - double t = eltrans->Weight(); - // Calculate r[i] = ve_j^T * elvect - for (int k = 0; k < elvect_size; k++) - { - coef[k + (i * elvect_size) + (j * rw.size() * elvect_size)] = vj_e[k] * rw[i] * t; - } - - // Calculate DS and store - CalcInverse(eltrans->Jacobian(), Jrt); - fe->CalcDShape(ip, DSh); - Mult(DSh, Jrt, DS); - for (int ii = 0; ii < dof; ++ii) - { - for (int jj = 0; jj < dim; ++jj) - { - index = jj + ii * dim; - DS_coef[index + (i * dof*dim) + (j * rw.size() * dof*dim)] = DS.Elem(ii,jj); - } - } - } - } -} + CAROM::Matrix const &V_v, const int rank, Vector &coef, Vector &DS_coef); void HyperelasticNLFIntegrator_ComputeReducedEQP(ParFiniteElementSpace *fesR, std::vector const &rw, std::vector const &qp, const IntegrationRule *ir, NeoHookeanModel *model, const Vector *x0, - CAROM::Matrix const &V_x, CAROM::Matrix const &V_v, CAROM::Vector const &x, const int rank, Vector &res) -{ - const int rxdim = V_x.numColumns(); - const int rvdim = V_v.numColumns(); - const int fomdim = V_v.numRows(); - MFEM_VERIFY(rw.size() == qp.size(), ""); - MFEM_VERIFY(x.dim() == rxdim, ""); - MFEM_VERIFY(V_x.numRows() == fesR->GetTrueVSize(), ""); - - MFEM_VERIFY(rank == 0, - "TODO: generalize to parallel. This uses full dofs in V, which has true dofs"); - - const int nqe = ir->GetWeights().Size(); - - ElementTransformation *eltrans; - DofTransformation *doftrans; - const FiniteElement *fe = NULL; - Array vdofs; - - res.SetSize(rxdim); - res = 0.0; - - int eprev = -1; - int dof = 0; - int dim = 0; - - // Get prolongation matrix - const Operator *P = fesR->GetProlongationMatrix(); - - // Vectors to be prolongated - CAROM::Vector *Vx_librom_temp = new CAROM::Vector(fomdim, true); - Vector *Vx_temp = new Vector(&((*Vx_librom_temp)(0)), fomdim); - Vector Vx(fomdim); - Vector vj(fomdim); - - // Prolongated vectors - Vector p_Vx(P->Height()); - Vector p_vj(P->Height()); - - // Element vectors - Vector Vx_e; - Vector vj_e; - - // Lift x, add x0 and prolongate result - V_x.mult(x, Vx_librom_temp); - add(*Vx_temp, *x0, Vx); - P->Mult(Vx, p_Vx); - - // Initialize nonlinear operator storage - // Assuming all elements have the same dim and n_dof - fe = fesR->GetFE(0); - dof = fe->GetDof(); - dim = fe->GetDim(); - DenseMatrix DSh(dof, dim); - DenseMatrix DS(dof, dim); - DenseMatrix Jrt(dim); - DenseMatrix Jpt(dim); - DenseMatrix P_f(dim); - DenseMatrix PMatI; // Extract element dofs - DenseMatrix PMatO; - Vector elvect(dof * dim); - PMatO.UseExternalData(elvect.GetData(), dof, dim); - - // For every basis vector - for (int j = 0; j < rvdim; ++j) - - { - // Get basis vector and prolongate - for (int k = 0; k < V_v.numRows(); ++k) - vj[k] = V_v(k, j); - P->Mult(vj, p_vj); - res[j] = 0.0; - - eprev = -1; - - // For every quadrature weight - for (int i = 0; i < rw.size(); ++i) // NOTE: i < 9 - { - const int e = qp[i] / nqe; // Element index - // Local (element) index of the quadrature point - const int qpi = qp[i] - (e * nqe); - const IntegrationPoint &ip = ir->IntPoint(qpi); - - if (e != eprev) // Update element transformation - { - doftrans = fesR->GetElementVDofs(e, vdofs); - fe = fesR->GetFE(e); - eltrans = fesR->GetElementTransformation(e); - - dof = fe->GetDof(); // Get number of dofs in element - dim = fe->GetDim(); - - if (doftrans) - { - MFEM_ABORT("TODO"); - } - - // Get element vectors - p_Vx.GetSubVector(vdofs, Vx_e); - p_vj.GetSubVector(vdofs, vj_e); - eprev = e; - } - - // Integration at ip - PMatI.UseExternalData(Vx_e.GetData(), dof, dim); - - elvect = 0.0; - - // Set integration point in the element transformation - eltrans->SetIntPoint(&ip); - model->SetTransformation(*eltrans); - - // Get the transformation weight - double t = eltrans->Weight(); - - // Compute action of nonlinear operator - CalcInverse(eltrans->Jacobian(), Jrt); - fe->CalcDShape(ip, DSh); - Mult(DSh, Jrt, DS); - MultAtB(PMatI, DS, Jpt); - model->EvalP(Jpt, P_f); - P_f *= (t * rw[i]); // NB: Not by ip.weight - AddMultABt(DS, P_f, PMatO); - - // Calculate r[i] = ve_j^T * elvect - for (int k = 0; k < elvect.Size(); k++) - { - res[j] += vj_e[k] * elvect[k]; - } - } - } -} + CAROM::Matrix const &V_x, CAROM::Matrix const &V_v, CAROM::Vector const &x, const int rank, Vector &res); void HyperelasticNLFIntegrator_ComputeReducedEQP_Fast(ParFiniteElementSpace *fesR, std::vector const &rw, std::vector const &qp, const IntegrationRule *ir, NeoHookeanModel *model, const Vector *x0, CAROM::Matrix const &V_x, CAROM::Matrix const &V_v, CAROM::Vector const &x, - Vector const &coef, Vector const &DS_coef, const int rank, Vector &res) -{ - - const int rxdim = V_x.numColumns(); - const int rvdim = V_v.numColumns(); - const int fomdim = V_x.numRows(); - MFEM_VERIFY(x.dim() == rxdim, ""); - MFEM_VERIFY(V_x.numRows() == fesR->GetTrueVSize(), ""); - - MFEM_VERIFY(rank == 0, - "TODO: generalize to parallel. This uses full dofs in V, which has true dofs"); - - const int nqe = ir->GetWeights().Size(); - - ElementTransformation *eltrans; - DofTransformation *doftrans; - const FiniteElement *fe = NULL; - Array vdofs; - - res.SetSize(rxdim); - res = 0.0; - - int eprev = -1; - int dof = 0; - int dim = 0; - - // Get prolongation matrix - const Operator *P = fesR->GetProlongationMatrix(); - - // Vectors to be prolongated - CAROM::Vector *Vx_librom_temp = new CAROM::Vector(fomdim, true); - Vector *Vx_temp = new Vector(&((*Vx_librom_temp)(0)), fomdim); - Vector Vx(fomdim); - - // Prolongated vectors - Vector p_Vx(P->Height()); - - // Element vectors - Vector Vx_e; - - // Lift x, add x0 and prolongate result - V_x.mult(x, Vx_librom_temp); - add(*Vx_temp, *x0, Vx); - P->Mult(Vx, p_Vx); - - // Initialize nonlinear operator storage - // Assuming all elements have the same dim and n_dof - fe = fesR->GetFE(0); - dof = fe->GetDof(); - dim = fe->GetDim(); - int index = 0; - DenseMatrix DS(dof, dim); - DenseMatrix Jpt(dim); - DenseMatrix P_f(dim); - DenseMatrix PMatI; // Extract element dofs - DenseMatrix PMatO; - Vector elvect(dof * dim); - PMatO.UseExternalData(elvect.GetData(), dof, dim); - - // For every basis vector - for (int j = 0; j < rvdim; ++j) - - { - - eprev = -1; - - // For every quadrature weight - for (int i = 0; i < qp.size(); ++i) // NOTE: i < 9 - { - const int e = qp[i] / nqe; // Element index - // Local (element) index of the quadrature point - const int qpi = qp[i] - (e * nqe); - const IntegrationPoint &ip = ir->IntPoint(qpi); - - if (e != eprev) // Update element transformation - { - doftrans = fesR->GetElementVDofs(e, vdofs); - fe = fesR->GetFE(e); - eltrans = fesR->GetElementTransformation(e); - - dof = fe->GetDof(); // Get number of dofs in element - dim = fe->GetDim(); - - if (doftrans) - { - MFEM_ABORT("TODO"); - } - - // Get element vectors - p_Vx.GetSubVector(vdofs, Vx_e); - eprev = e; - } + Vector const &coef, Vector const &DS_coef, const int rank, Vector &res); - // Integration at ip - elvect = 0.0; - PMatI.UseExternalData(Vx_e.GetData(), dof, dim); +bool fileExists(const std::string &filename); - // Set integration point in the element transformation - eltrans->SetIntPoint(&ip); - model->SetTransformation(*eltrans); - - for (int ii = 0; ii < dof; ++ii) - { - for (int jj = 0; jj < dim; ++jj) - { - index = jj + ii * dim; - DS.Elem(ii,jj) = DS_coef[index + (i * elvect.Size()) + (j * qp.size() * elvect.Size())]; - } - } - - MultAtB(PMatI, DS, Jpt); - model->EvalP(Jpt, P_f); - AddMultABt(DS, P_f, PMatO); - - // Calculate r[i] = ve_j^T * elvect - // coef is size len(vdofs) * rvdim * rw.size - for (int k = 0; k < elvect.Size(); k++) - { - res[j] += coef[k + (i * elvect.Size()) + (j * qp.size() * elvect.Size())] * elvect[k]; - } - } - } -} - -bool fileExists(const std::string &filename) -{ - std::ifstream file(filename); - return file.good(); -} - -void save_CAROM_Vector(const CAROM::Vector &vector, const std::string &filename) -{ - std::ofstream file(filename); - if (file.is_open()) - { - for (int i = 0; i < vector.dim(); ++i) - { - file << vector(i) << "\n"; - } - file.close(); - std::cout << "Vector saved as " << filename << " successfully." << std::endl; - } - else - { - std::cerr << "Error: Unable to open file " << filename << " for writing." << std::endl; - } -} +void save_CAROM_Vector(const CAROM::Vector &vector, const std::string &filename); void ComputeElementRowOfG(const IntegrationRule *ir, Array const &vdofs, Vector const &ve_j, NeoHookeanModel *model, Vector const &elfun, - FiniteElement const &fe, ElementTransformation &Trans, Vector &r) -{ - MFEM_VERIFY(r.Size() == ir->GetNPoints(), ""); - int dof = fe.GetDof(); // Get number of dofs in element - int dim = fe.GetDim(); - - // Initialize nonlinear operator matrices (there is probably a better way) - DenseMatrix DSh(dof, dim); - DenseMatrix DS(dof, dim); - DenseMatrix Jrt(dim); - DenseMatrix Jpt(dim); - DenseMatrix P(dim); - DenseMatrix PMatI; // Extract element dofs - PMatI.UseExternalData(elfun.GetData(), dof, dim); - DenseMatrix PMatO; - Vector elvect(dof * dim); - PMatO.UseExternalData(elvect.GetData(), dof, dim); - - model->SetTransformation(Trans); - - // For each integration point - for (int i = 0; i < ir->GetNPoints(); i++) - { - elvect = 0.0; - // Get integration point - const IntegrationPoint &ip = ir->IntPoint(i); - - // Set integration point in the element transformation - Trans.SetIntPoint(&ip); - - // Get the transformation weight - double t = Trans.Weight(); - - // Compute action of nonlinear operator - CalcInverse(Trans.Jacobian(), Jrt); - fe.CalcDShape(ip, DSh); - Mult(DSh, Jrt, DS); - MultAtB(PMatI, DS, Jpt); - model->EvalP(Jpt, P); - P *= t; // NB: Not by ip.weight - AddMultABt(DS, P, PMatO); - - r[i] = 0.0; - - // Calculate r[i] = ve_j^T * elvect - for (int k = 0; k < elvect.Size(); k++) - { - r[i] += ve_j[k] * elvect[k]; - } - } -} + FiniteElement const &fe, ElementTransformation &Trans, Vector &r); void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, CAROM::Vector const &w, CAROM::Matrix &Gt, - CAROM::Vector &sol) -{ - 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); - - CAROM::Vector rhs_lb(rhs_ub); - CAROM::Vector rhs_Gw(rhs_ub); - - const double delta = 1.0e-11; - for (int i = 0; i < rhs_ub.dim(); ++i) - { - rhs_lb(i) -= delta; - rhs_ub(i) += delta; - } - - // nnls.normalize_constraints(Gt, rhs_lb, rhs_ub); - nnls.solve_parallel_with_scalapack(Gt, rhs_lb, rhs_ub, sol); - - int nnz = 0; - for (int i = 0; i < sol.dim(); ++i) - { - if (sol(i) != 0.0) - { - nnz++; - } - } - - cout << rank << ": Number of nonzeros in NNLS solution: " << nnz - << ", out of " << sol.dim() << endl; - - MPI_Allreduce(MPI_IN_PLACE, &nnz, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); - - if (rank == 0) - cout << "Global number of nonzeros in NNLS solution: " << nnz << endl; - - // Check residual of NNLS solution - CAROM::Vector res(Gt.numColumns(), false); - Gt.transposeMult(sol, res); - - const double normGsol = res.norm(); - const double normRHS = rhs_Gw.norm(); - - res -= rhs_Gw; - const double relNorm = res.norm() / std::max(normGsol, normRHS); - cout << rank << ": relative residual norm for NNLS solution of Gs = Gw: " << relNorm << endl; -} + CAROM::Vector &sol); // Compute EQP solution from constraints on snapshots. void SetupEQP_snapshots(const IntegrationRule *ir0, const int rank, @@ -568,246 +91,16 @@ void SetupEQP_snapshots(const IntegrationRule *ir0, const int rank, const CAROM::Matrix *BX_snapshots, const bool precondition, const double nnls_tol, const int maxNNLSnnz, NeoHookeanModel *model, - CAROM::Vector &sol, CAROM::Vector *window_ids) -{ - const int nqe = ir0->GetNPoints(); - const int ne = fespace_X->GetNE(); - const int NB = BV->numColumns(); - const int NQ = ne * nqe; - const int nsnap = BX_snapshots->numColumns(); - - MFEM_VERIFY(nsnap == BX_snapshots->numColumns() || - nsnap + nsets == BX_snapshots->numColumns(), // Q: nsets? - ""); - MFEM_VERIFY(BV->numRows() == BX_snapshots->numRows(), ""); - MFEM_VERIFY(BV->numRows() == fespace_X->GetTrueVSize(), ""); - - const bool skipFirstW = (nsnap + nsets == BX_snapshots->numColumns()); - - // Compute G of size (NB * nsnap) x NQ, but only store its transpose Gt. - CAROM::Matrix Gt(NQ, NB * nsnap, true); - int current_size = 0; - int i_start = 0; - int outer_loop_length = 0; - if (!window_ids) - { - current_size = nsnap; - outer_loop_length = 1; - } - else - { - outer_loop_length = window_ids->dim() - 1; - } - - // For 0 <= j < NB, 0 <= i < nsnap, 0 <= e < ne, 0 <= m < nqe, - // G(j + (i*NB), (e*nqe) + m) - // is the coefficient of v_j^T M(p_i) V 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 x_i(BX_snapshots->numRows()); - Vector v_j(BV->numRows()); - - Vector r(nqe); - - int skip = 0; - const int nsnapPerSet = nsnap / nsets; - if (skipFirstW) - { - MFEM_VERIFY(nsets * nsnapPerSet == nsnap, ""); - skip = 1; - } - - // Get prolongation matrix - const Operator *P = fespace_X->GetProlongationMatrix(); - if (!P) - { - MFEM_ABORT("P is null, generalize to serial case") - } - - // Outer loop for time windowing - for (int oi = 0; oi < outer_loop_length; ++oi) - { - if (window_ids) - { - i_start = window_ids->item(oi); - current_size = window_ids->item(oi + 1) - window_ids->item(oi) + 1; - cout << "i_start = " << i_start << endl; - Gt.setSize(NQ, NB * current_size); - } - - // For every snapshot in batch - for (int i = i_start; i < current_size; ++i) - { - // Set the sampled dofs from the snapshot matrix - for (int j = 0; j < BX_snapshots->numRows(); ++j) - x_i[j] = (*BX_snapshots)(j, i); - - // Get prolongated dofs - Vector px_i; - Vector elfun; - - px_i.SetSize(P->Height()); - P->Mult(x_i, px_i); - - if (skipFirstW && i > 0 && i % nsnapPerSet == 0) - skip++; - - // For each basis vector - for (int j = 0; j < NB; ++j) - { - - // Get basis vector - for (int k = 0; k < BV->numRows(); ++k) - v_j[k] = (*BV)(k, j); - - // Get prolongated dofs - Vector pv_j; - Vector ve_j; - - pv_j.SetSize(P->Height()); - P->Mult(v_j, pv_j); - - // TODO: is it better to make the element loop the outer loop? - // For each element - for (int e = 0; e < ne; ++e) - { - // Get element and its dofs and transformation. - Array vdofs; - DofTransformation *doftrans = fespace_X->GetElementVDofs(e, vdofs); - const FiniteElement &fe = *fespace_X->GetFE(e); - ElementTransformation *eltrans = fespace_X->GetElementTransformation(e); - px_i.GetSubVector(vdofs, elfun); - pv_j.GetSubVector(vdofs, ve_j); - if (doftrans) - { - MFEM_ABORT("Doftrans is true, make corresponding edits") - } - - // Compute the row of G corresponding to element e, store in r - ComputeElementRowOfG(ir0, vdofs, ve_j, model, elfun, fe, *eltrans, r); - - for (int m = 0; m < nqe; ++m) - Gt((e * nqe) + m, j + (i * NB)) = r[m]; - } - } - - if (precondition) - { - MFEM_ABORT("TODO"); - } - } // Loop (i) over snapshots - - Array const &w_el = ir0->GetWeights(); - MFEM_VERIFY(w_el.Size() == nqe, ""); - - CAROM::Vector w(ne * nqe, true); - - for (int i = 0; i < ne; ++i) - { - for (int j = 0; j < nqe; ++j) - w((i * nqe) + j) = w_el[j]; - } - - SolveNNLS(rank, nnls_tol, maxNNLSnnz, w, Gt, sol); - if (window_ids) - { - save_CAROM_Vector(sol, "sol_" + std::to_string(oi) + ".csv"); - } - } -} + CAROM::Vector &sol, CAROM::Vector *window_ids); void WriteMeshEQP(ParMesh *pmesh, const int myid, const int nqe, - CAROM::Vector const &eqpSol) -{ - // Find the elements with quadrature points in eqpSol. - std::set elements; - Vector elemCount(pmesh->GetNE()); - elemCount = 0.0; - - for (int i = 0; i < eqpSol.dim(); ++i) - { - if (eqpSol(i) > 1.0e-12) - { - const int e = i / nqe; // Element index - elements.insert(e); - elemCount[e] += 1.0; - } - } - - // Empty sets, since EQP on samples inside elements. - std::set faces, edges, vertices; - CAROM::SampleVisualization(pmesh, elements, elements, faces, edges, - vertices, "EQPvis", &elemCount); -} + CAROM::Vector const &eqpSol); // Function to compute the indices at which to switch time windows -void get_window_ids(int n_step, int n_window, CAROM::Vector *ids) -{ - int window_size = std::round(n_step / n_window); - - int res = n_step - (n_window * (window_size - 1) + 1); - int res_up = res - n_window; - int ctr = -1; - - for (int i = 1; i <= n_window + 1; ++i) - { - ids->item(i - 1) = (i - 1) * window_size + 1 - (i - 1); - if (res_up > 0) - { - if (i <= res - res_up) - { - ctr += 1; - } - if (i > 2 && i <= res_up + 1) - { - ctr += 1; - } - } - else - { - if (i <= res + 1) - { - ctr += 1; - } - } - ids->item(i - 1) += ctr; - } -} - -void load_CAROM_vector(const std::string &filename, CAROM::Vector &vector) -{ - std::ifstream file(filename); - std::string line; - std::vector data; - - while (std::getline(file, line)) - { - std::stringstream ss(line); - std::string value; +void get_window_ids(int n_step, int n_window, CAROM::Vector *ids); - while (std::getline(ss, value, ',')) - { - data.push_back(std::stod(value)); - } - } +void load_CAROM_vector(const std::string &filename, CAROM::Vector &vector); - // Set the size of the vector - int size = data.size(); - vector.setSize(size); +void get_EQPsol(const int current_window, CAROM::Vector *load_eqpsol); - // Copy the data into the vector - for (int i = 0; i < size; ++i) - { - vector(i) = data[i]; - } -} - -void get_EQPsol(const int current_window, CAROM::Vector *load_eqpsol) -{ - string filename = "sol_" + std::to_string(current_window) + ".csv"; - if (fileExists(filename)) - { - load_CAROM_vector(filename, *load_eqpsol); - } -} +#endif