diff --git a/.gitignore b/.gitignore index bbd2846d02..f87ea59b9a 100644 --- a/.gitignore +++ b/.gitignore @@ -296,6 +296,7 @@ miniapps/tools/display-basis miniapps/tools/load-dc miniapps/tools/convert-dc miniapps/tools/lor-transfer +miniapps/tools/plor-transfer miniapps/tools/get-values miniapps/tools/check-tmop-metric miniapps/tools/tmop-metric-magnitude diff --git a/CHANGELOG b/CHANGELOG index 30ec22e9d2..f7011929a4 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -45,6 +45,10 @@ New and updated examples and miniapps - Moved the distance solver methods from miniapps/shifted to miniapps/common. +- Added a new parallel LOR transfer miniapp, miniapps/tools/plor-transfer, which + mirrors the functionality of the serial LOR transfer miniapp, + miniapps/tools/lor-transfer + Meshing improvements -------------------- - Added support for free connectivity of NURBS patches allowing for more complex @@ -81,6 +85,9 @@ Discretization improvements - Added support for p-refined meshes in FindPointsGSLIB. +- Support for parallel transfer of H1 fields using the low-order refined (LOR) + transfer operators in L2ProjectionGridTransfer + Linear and nonlinear solvers ---------------------------- - Updated interface to MUMPS direct solver to support multiple right-hand diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 68dab72d5d..9ca9e5f3cc 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -546,6 +546,234 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_) : L2Projection(fes_ho_, fes_lor_) { + std::unique_ptr R_mat, M_LH_mat; + std::tie(R_mat, M_LH_mat) = ComputeSparseRAndM_LH(); + + FiniteElementSpace fes_ho_scalar(fes_ho.GetMesh(), fes_ho.FEColl(), 1); + FiniteElementSpace fes_lor_scalar(fes_lor.GetMesh(), fes_lor.FEColl(), 1); + + const SparseMatrix *P_ho = fes_ho_scalar.GetConformingProlongation(); + const SparseMatrix *P_lor = fes_lor_scalar.GetConformingProlongation(); + + if (P_ho || P_lor) + { + if (P_ho && P_lor) + { + R_mat.reset(RAP(*P_lor, *R_mat, *P_ho)); + M_LH_mat.reset(RAP(*P_lor, *M_LH_mat, *P_ho)); + } + else if (P_ho) + { + R_mat.reset(mfem::Mult(*R_mat, *P_ho)); + M_LH_mat.reset(mfem::Mult(*M_LH_mat, *P_ho)); + } + else // P_lor != nullptr + { + R_mat.reset(mfem::Mult(*P_lor, *R_mat)); + M_LH_mat.reset(mfem::Mult(*P_lor, *M_LH_mat)); + } + } + + SparseMatrix *RTxM_LH_mat = TransposeMult(*R_mat, *M_LH_mat); + precon.reset(new DSmoother(*RTxM_LH_mat)); + + // Set ownership + RTxM_LH.reset(RTxM_LH_mat); + R = std::move(R_mat); + M_LH = std::move(M_LH_mat); + + SetupPCG(); +} + +#ifdef MFEM_USE_MPI + +L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( + const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor) + : L2Projection(pfes_ho, pfes_lor), + pcg(pfes_ho.GetComm()) +{ + std::tie(R, M_LH) = ComputeSparseRAndM_LH(); + + ParFiniteElementSpace pfes_ho_scalar(pfes_ho.GetParMesh(), + pfes_ho.FEColl(), 1); + ParFiniteElementSpace pfes_lor_scalar(pfes_lor.GetParMesh(), + pfes_lor.FEColl(), 1); + + HypreParMatrix R_local = HypreParMatrix(pfes_ho.GetComm(), + pfes_lor_scalar.GlobalVSize(), + pfes_ho_scalar.GlobalVSize(), + pfes_lor_scalar.GetDofOffsets(), + pfes_ho_scalar.GetDofOffsets(), + static_cast(R.get())); + HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho.GetComm(), + pfes_lor_scalar.GlobalVSize(), + pfes_ho_scalar.GlobalVSize(), + pfes_lor_scalar.GetDofOffsets(), + pfes_ho_scalar.GetDofOffsets(), + static_cast(M_LH.get())); + + HypreParMatrix *R_mat = RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), + &R_local, pfes_ho_scalar.Dof_TrueDof_Matrix()); + HypreParMatrix *M_LH_mat = RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), + &M_LH_local, pfes_ho_scalar.Dof_TrueDof_Matrix()); + + std::unique_ptr R_T(R_mat->Transpose()); + HypreParMatrix *RTxM_LH_mat = ParMult(R_T.get(), M_LH_mat, true); + + HypreBoomerAMG *amg = new HypreBoomerAMG(*RTxM_LH_mat); + amg->SetPrintLevel(0); + + R.reset(R_mat); + M_LH.reset(M_LH_mat); + RTxM_LH.reset(RTxM_LH_mat); + precon.reset(amg); + + SetupPCG(); + pcg.SetPreconditioner(*precon); + pcg.SetOperator(*RTxM_LH); +} + +#endif + +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetupPCG() +{ + // Basic PCG solver setup + pcg.SetPrintLevel(0); + // pcg.SetPrintLevel(IterativeSolver::PrintLevel().Summary()); + pcg.SetMaxIter(1000); + // initial values for relative and absolute tolerance + pcg.SetRelTol(1e-13); + pcg.SetAbsTol(1e-13); + pcg.SetPreconditioner(*precon); + pcg.SetOperator(*RTxM_LH); +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( + const Vector& x, Vector& y) const +{ + Vector X(fes_ho.GetTrueVSize()); + Vector X_dim(R->Width()); + + Vector Y_dim(R->Height()); + Vector Y(fes_lor.GetTrueVSize()); + + Array vdofs_list; + + GetTDofs(fes_ho, x, X); + + for (int d = 0; d < fes_ho.GetVDim(); ++d) + { + TDofsListByVDim(fes_ho, d, vdofs_list); + X.GetSubVector(vdofs_list, X_dim); + R->Mult(X_dim, Y_dim); + TDofsListByVDim(fes_lor, d, vdofs_list); + Y.SetSubVector(vdofs_list, Y_dim); + } + + SetFromTDofs(fes_lor, Y, y); +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( + const Vector& x, Vector& y) const +{ + Vector X(fes_lor.GetTrueVSize()); + Vector X_dim(R->Height()); + + Vector Y_dim(R->Width()); + Vector Y(fes_ho.GetTrueVSize()); + + Array vdofs_list; + + GetTDofsTranspose(fes_lor, x, X); + + for (int d = 0; d < fes_ho.GetVDim(); ++d) + { + TDofsListByVDim(fes_lor, d, vdofs_list); + X.GetSubVector(vdofs_list, X_dim); + R->MultTranspose(X_dim, Y_dim); + TDofsListByVDim(fes_ho, d, vdofs_list); + Y.SetSubVector(vdofs_list, Y_dim); + } + + SetFromTDofsTranspose(fes_ho, Y, y); +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( + const Vector& x, Vector& y) const +{ + Vector X(fes_lor.GetTrueVSize()); + Vector X_dim(M_LH->Height()); + Vector Xbar(pcg.Width()); + + Vector Y_dim(pcg.Height()); + Vector Y(fes_ho.GetTrueVSize()); + + Array vdofs_list; + + GetTDofs(fes_lor, x, X); + + for (int d = 0; d < fes_ho.GetVDim(); ++d) + { + TDofsListByVDim(fes_lor, d, vdofs_list); + X.GetSubVector(vdofs_list, X_dim); + // Compute y = P x = (R^T M_LH)^(-1) M_LH^T X = (R^T M_LH)^(-1) Xbar + M_LH->MultTranspose(X_dim, Xbar); + Y_dim = 0.0; + pcg.Mult(Xbar, Y_dim); + TDofsListByVDim(fes_ho, d, vdofs_list); + Y.SetSubVector(vdofs_list, Y_dim); + } + + SetFromTDofs(fes_ho, Y, y); +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( + const Vector& x, Vector& y) const +{ + Vector X(fes_ho.GetTrueVSize()); + Vector X_dim(pcg.Width()); + Vector Xbar(pcg.Height()); + + Vector Y_dim(M_LH->Height()); + Vector Y(fes_lor.GetTrueVSize()); + + Array vdofs_list; + + GetTDofsTranspose(fes_ho, x, X); + + for (int d = 0; d < fes_ho.GetVDim(); ++d) + { + TDofsListByVDim(fes_ho, d, vdofs_list); + X.GetSubVector(vdofs_list, X_dim); + // Compute y = P^T x = M_LH (R^T M_LH)^(-1) X = M_LH Xbar + Xbar = 0.0; + pcg.Mult(X_dim, Xbar); + M_LH->Mult(Xbar, Y_dim); + TDofsListByVDim(fes_lor, d, vdofs_list); + Y.SetSubVector(vdofs_list, Y_dim); + } + + SetFromTDofsTranspose(fes_lor, Y, y); +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_) +{ + pcg.SetRelTol(p_rtol_); +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_) +{ + pcg.SetAbsTol(p_atol_); +} + +std::pair< +std::unique_ptr, +std::unique_ptr> + L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() +{ + std::pair, + std::unique_ptr> r_and_mlh; + Mesh* mesh_ho = fes_ho.GetMesh(); Mesh* mesh_lor = fes_lor.GetMesh(); int nel_ho = mesh_ho->GetNE(); @@ -553,7 +781,7 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( int ndof_lor = fes_lor.GetNDofs(); // If the local mesh is empty, skip all computations - if (nel_ho == 0) { return; } + if (nel_ho == 0) { return {nullptr, nullptr}; } const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms(); @@ -611,18 +839,26 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( } } // DOF by DOF inverse of non-zero entries - for (int i = 0; i < ndof_lor; ++i) - { - ML_inv[i] = 1.0 / ML_inv[i]; - } + LumpedMassInverse(ML_inv); // Compute sparsity pattern for R = M_L^(-1) M_LH and allocate - AllocR(); + r_and_mlh.first = AllocR(); // Allocate M_LH (same sparsity pattern as R) // L refers to the low-order refined mesh (DOFs correspond to rows) // H refers to the higher-order mesh (DOFs correspond to columns) - M_LH = SparseMatrix(R.GetI(), R.GetJ(), NULL, - R.Height(), R.Width(), false, true, true); + Memory I(r_and_mlh.first->Height() + 1); + for (int icol = 0; icol < r_and_mlh.first->Height() + 1; ++icol) + { + I[icol] = r_and_mlh.first->GetI()[icol]; + } + Memory J(r_and_mlh.first->NumNonZeroElems()); + for (int jcol = 0; jcol < r_and_mlh.first->NumNonZeroElems(); ++jcol) + { + J[jcol] = r_and_mlh.first->GetJ()[jcol]; + } + r_and_mlh.second = std::unique_ptr(new SparseMatrix( + I, J, NULL, + r_and_mlh.first->Height(), r_and_mlh.first->Width(), true, true, true)); IntegrationPointTransformation ip_tr; IsoparametricTransformation& emb_tr = ip_tr.Transf; @@ -667,131 +903,118 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( } Array dofs_ho(nedof_ho); fes_ho.GetElementDofs(iho, dofs_ho); - M_LH.AddSubMatrix(dofs_lor, dofs_ho, M_LH_el); - R.AddSubMatrix(dofs_lor, dofs_ho, R_el); + r_and_mlh.second->AddSubMatrix(dofs_lor, dofs_ho, M_LH_el); + r_and_mlh.first->AddSubMatrix(dofs_lor, dofs_ho, R_el); } } - // Create PCG solver - RTxM_LH = TransposeMult(R, M_LH); - pcg.SetPrintLevel(0); - pcg.SetMaxIter(1000); - // initial values for relative and absolute tolerance - SetRelTol(1e-13); - SetAbsTol(1e-13); - Ds = DSmoother(*RTxM_LH); - pcg.SetPreconditioner(Ds); - pcg.SetOperator(*RTxM_LH); + return r_and_mlh; } -L2ProjectionGridTransfer::L2ProjectionH1Space::~L2ProjectionH1Space() +void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofs( + const FiniteElementSpace& fes, const Vector& x, Vector& X) const { - delete RTxM_LH; -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( - const Vector& x, Vector& y) const -{ - int vdim = fes_ho.GetVDim(); - const int ndof_ho = fes_ho.GetNDofs(); - const int ndof_lor = fes_lor.GetNDofs(); - Array dofs_ho(ndof_ho); - Array dofs_lor(ndof_lor); - Vector x_dim(ndof_ho); - Vector y_dim(ndof_lor); - - for (int d = 0; d < vdim; ++d) + const Operator* res = fes.GetRestrictionOperator(); + if (res) { - fes_ho.GetVDofs(d, dofs_ho); - fes_lor.GetVDofs(d, dofs_lor); - x.GetSubVector(dofs_ho, x_dim); - R.Mult(x_dim, y_dim); - y.SetSubVector(dofs_lor, y_dim); + res->Mult(x, X); + } + else + { + X = x; } } -void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofs( + const FiniteElementSpace& fes, const Vector &X, Vector& x) const { - int vdim = fes_ho.GetVDim(); - const int ndof_ho = fes_ho.GetNDofs(); - const int ndof_lor = fes_lor.GetNDofs(); - Array dofs_ho(ndof_ho); - Array dofs_lor(ndof_lor); - Vector x_dim(ndof_lor); - Vector y_dim(ndof_ho); - - for (int d = 0; d < vdim; ++d) + const Operator* P = fes.GetProlongationMatrix(); + if (P) { - fes_ho.GetVDofs(d, dofs_ho); - fes_lor.GetVDofs(d, dofs_lor); - x.GetSubVector(dofs_lor, x_dim); - R.MultTranspose(x_dim, y_dim); - y.SetSubVector(dofs_ho, y_dim); + P->Mult(X, x); + } + else + { + x = X; } } -void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofsTranspose( + const FiniteElementSpace& fes, const Vector& x, Vector& X) const { - int vdim = fes_ho.GetVDim(); - const int ndof_ho = fes_ho.GetNDofs(); - const int ndof_lor = fes_lor.GetNDofs(); - Array dofs_ho(ndof_ho); - Array dofs_lor(ndof_lor); - Vector x_dim(ndof_lor); - Vector y_dim(ndof_ho); - Vector xbar(ndof_ho); - - for (int d = 0; d < vdim; ++d) + const Operator* P = fes.GetProlongationMatrix(); + if (P) { - fes_lor.GetVDofs(d, dofs_lor); - x.GetSubVector(dofs_lor, x_dim); - // Compute y = P x = (R^T M_LH)^(-1) M_LH^T x = (R^T M_LH)^(-1) xbar - M_LH.MultTranspose(x_dim, xbar); - y_dim = 0.0; - pcg.Mult(xbar, y_dim); - fes_ho.GetVDofs(d, dofs_ho); - y.SetSubVector(dofs_ho, y_dim); + P->MultTranspose(x, X); + } + else + { + X = x; } } -void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofsTranspose( + const FiniteElementSpace& fes, const Vector &X, Vector& x) const { - int vdim = fes_ho.GetVDim(); - const int ndof_ho = fes_ho.GetNDofs(); - const int ndof_lor = fes_lor.GetNDofs(); - Array dofs_ho(ndof_ho); - Array dofs_lor(ndof_lor); - Vector x_dim(ndof_ho); - Vector y_dim(ndof_lor); - Vector xbar(ndof_ho); - - for (int d = 0; d < vdim; ++d) + const Operator *R_op = fes.GetRestrictionOperator(); + if (R_op) { - fes_ho.GetVDofs(d, dofs_ho); - x.GetSubVector(dofs_ho, x_dim); - // Compute y = P^T x = M_LH (R^T M_LH)^(-1) x = M_LH xbar - xbar = 0.0; - pcg.Mult(x_dim, xbar); - M_LH.Mult(xbar, y_dim); - fes_lor.GetVDofs(d, dofs_lor); - y.SetSubVector(dofs_lor, y_dim); + R_op->MultTranspose(X, x); + } + else + { + x = X; } } -void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_) +void L2ProjectionGridTransfer::L2ProjectionH1Space::TDofsListByVDim( + const FiniteElementSpace& fes, int vdim, Array& vdofs_list) const { - pcg.SetRelTol(p_rtol_); + const SparseMatrix *R_mat = fes.GetRestrictionMatrix(); + if (R_mat) + { + Array x_vdofs_list(fes.GetNDofs()); + Array x_vdofs_marker(fes.GetVSize()); + Array X_vdofs_marker(fes.GetTrueVSize()); + fes.GetVDofs(vdim, x_vdofs_list); + FiniteElementSpace::ListToMarker(x_vdofs_list, fes.GetVSize(), x_vdofs_marker); + R_mat->BooleanMult(x_vdofs_marker, X_vdofs_marker); + FiniteElementSpace::MarkerToList(X_vdofs_marker, vdofs_list); + } + else + { + vdofs_list.SetSize(fes.GetNDofs()); + fes.GetVDofs(vdim, vdofs_list); + } } -void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_) +void L2ProjectionGridTransfer::L2ProjectionH1Space::LumpedMassInverse( + Vector& ML_inv) const { - pcg.SetAbsTol(p_atol_); + Vector ML_inv_full(fes_lor.GetVSize()); + // set ML_inv on dofs for vdim = 0 + Array vdofs_list(fes_lor.GetNDofs()); + fes_lor.GetVDofs(0, vdofs_list); + ML_inv_full.SetSubVector(vdofs_list, ML_inv); + + Vector ML_inv_true(fes_lor.GetTrueVSize()); + const Operator *P = fes_lor.GetProlongationMatrix(); + if (P) { P->MultTranspose(ML_inv_full, ML_inv_true); } + else { ML_inv_true = ML_inv_full; } + + for (int i = 0; i < ML_inv_true.Size(); ++i) + { + ML_inv_true[i] = 1.0 / ML_inv_true[i]; + } + + if (P) { P->Mult(ML_inv_true, ML_inv_full); } + else { ML_inv_full = ML_inv_true; } + + ML_inv_full.GetSubVector(vdofs_list, ML_inv); } -void L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() +std::unique_ptr +L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() { const Table& elem_dof_ho = fes_ho.GetElementToDofTable(); const Table& elem_dof_lor = fes_lor.GetElementToDofTable(); @@ -871,11 +1094,13 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() dof_lor_dof_ho.SortRows(); double* data = Memory(dof_dofI[ndof_lor]); - R = SparseMatrix(dof_dofI, dof_dofJ, data, ndof_lor, ndof_ho, - true, true, true); - R = 0.0; + std::unique_ptr R_local(new SparseMatrix( + dof_dofI, dof_dofJ, data, ndof_lor, ndof_ho, true, true, true)); + (*R_local) = 0.0; dof_lor_dof_ho.LoseData(); + + return R_local; } L2ProjectionGridTransfer::~L2ProjectionGridTransfer() @@ -905,7 +1130,20 @@ void L2ProjectionGridTransfer::BuildF() if (!force_l2_space && dom_fes.FEColl()->GetContType() == FiniteElementCollection::CONTINUOUS) { - F = new L2ProjectionH1Space(dom_fes, ran_fes); + if (!Parallel()) + { + F = new L2ProjectionH1Space(dom_fes, ran_fes); + } + else + { +#ifdef MFEM_USE_MPI + const mfem::ParFiniteElementSpace& dom_pfes = + static_cast(dom_fes); + const mfem::ParFiniteElementSpace& ran_pfes = + static_cast(ran_fes); + F = new L2ProjectionH1Space(dom_pfes, ran_pfes); +#endif + } } else { diff --git a/fem/transfer.hpp b/fem/transfer.hpp index a0a6988788..553f037cc6 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -180,9 +180,15 @@ protected: public: virtual void Prolongate(const Vector& x, Vector& y) const = 0; virtual void ProlongateTranspose(const Vector& x, Vector& y) const = 0; - /// Sets relative tolerance and absolute tolerance in preconditioned - /// conjugate gradient solver. Only used for H1 spaces. + /// @brief Sets relative tolerance in preconditioned conjugate gradient + /// solver. + /// + /// Only used for H1 spaces. virtual void SetRelTol(double p_rtol_) = 0; + /// @brief Sets absolute tolerance in preconditioned conjugate gradient + /// solver. + /// + /// Only used for H1 spaces. virtual void SetAbsTol(double p_atol_) = 0; protected: const FiniteElementSpace& fes_ho; @@ -249,29 +255,22 @@ protected: /// conservative left-inverse prolongation operation. This functionality /// is also provided as an Operator by L2Prolongation. virtual void ProlongateTranspose(const Vector& x, Vector& y) const; - virtual void SetRelTol(double p_rtol_) {} - virtual void SetAbsTol(double p_atol_) {} + virtual void SetRelTol(double p_rtol_) { } ///< No-op. + virtual void SetAbsTol(double p_atol_) { } ///< No-op. }; - /** Class for projection operator between a H1 high-order finite element - space on a coarse mesh, and a H1 low-order finite element space on a - refined mesh (LOR). */ + /** Projection operator between a H1 high-order finite element space on a + coarse mesh, and a H1 low-order finite element space on a refined mesh + (LOR). */ class L2ProjectionH1Space : public L2Projection { - // The restriction operator is represented as a SparseMatrix R. The - // prolongation operator is a dense matrix computed as the inverse of (R^T - // M_L R), and hence, is not stored. - SparseMatrix R; - // Used to compute P = (RTxM_LH)^(-1) M_LH^T - SparseMatrix M_LH; - SparseMatrix* RTxM_LH; - CGSolver pcg; - DSmoother Ds; - public: - L2ProjectionH1Space(const FiniteElementSpace& fes_ho_, - const FiniteElementSpace& fes_lor_); - virtual ~L2ProjectionH1Space(); + L2ProjectionH1Space(const FiniteElementSpace &fes_ho_, + const FiniteElementSpace &fes_lor_); +#ifdef MFEM_USE_MPI + L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_, + const ParFiniteElementSpace &pfes_lor_); +#endif /// Maps x, primal field coefficients defined on a coarse mesh /// with a higher order H1 finite element space, to y, primal /// field coefficients defined on a refined mesh with a low order H1 @@ -305,11 +304,51 @@ protected: virtual void ProlongateTranspose(const Vector& x, Vector& y) const; virtual void SetRelTol(double p_rtol_); virtual void SetAbsTol(double p_atol_); - private: - /// Computes sparsity pattern and initializes R matrix. Based on - /// BilinearForm::AllocMat() except maps between HO elements and LOR - /// elements. - void AllocR(); + protected: + /// Sets up the PCG solver (sets parameters, operator, and preconditioner) + void SetupPCG(); + /// Computes on-rank R and M_LH matrices. + std::pair, + std::unique_ptr> ComputeSparseRAndM_LH(); + /// @brief Recovers vector of tdofs given a vector of dofs and a finite + /// element space + void GetTDofs(const FiniteElementSpace& fes, const Vector& x, Vector& X) const; + /// Sets dof values given a vector of tdofs and a finite element space + void SetFromTDofs(const FiniteElementSpace& fes, + const Vector& X, + Vector& x) const; + /// @brief Recovers a vector of dual field coefficients on the tdofs given + /// a vector of dual coefficients and a finite element space + void GetTDofsTranspose(const FiniteElementSpace& fes, + const Vector& x, + Vector& X) const; + /// @brief Sets dual field coefficients given a vector of dual field + /// coefficients on the tdofs and a finite element space + void SetFromTDofsTranspose(const FiniteElementSpace& fes, + const Vector& X, + Vector& x) const; + /// @brief Fills the vdofs_list array with a list of vdofs for a given + /// vdim and a given finite element space + void TDofsListByVDim(const FiniteElementSpace& fes, + int vdim, + Array& vdofs_list) const; + /// Returns the inverse of an on-rank lumped mass matrix + void LumpedMassInverse(Vector& ML_inv) const; + /// @brief Computes sparsity pattern and initializes R matrix. + /// + /// Based on BilinearForm::AllocMat(), except maps between coarse HO + /// elements and refined LOR elements. + std::unique_ptr AllocR(); + + CGSolver pcg; + std::unique_ptr precon; + // The restriction operator is represented as an Operator R. The + // prolongation operator is a dense matrix computed as the inverse of (R^T + // M_L R), and hence, is not stored. + std::unique_ptr R; + // Used to compute P = (RT*M_LH)^(-1) M_LH^T + std::unique_ptr M_LH; + std::unique_ptr RTxM_LH; }; /** Mass-conservative prolongation operator going in the opposite direction diff --git a/miniapps/tools/CMakeLists.txt b/miniapps/tools/CMakeLists.txt index f746d2fd2c..cf6ba6966c 100644 --- a/miniapps/tools/CMakeLists.txt +++ b/miniapps/tools/CMakeLists.txt @@ -28,5 +28,10 @@ add_mfem_miniapp(convert-dc add_mfem_miniapp(lor-transfer MAIN lor-transfer.cpp LIBRARIES mfem) +if (MFEM_USE_MPI) + add_mfem_miniapp(plor-transfer + MAIN plor-transfer.cpp LIBRARIES mfem) +endif() + add_mfem_miniapp(check-tmop-metric MAIN check-tmop-metric.cpp LIBRARIES mfem) diff --git a/miniapps/tools/lor-transfer.cpp b/miniapps/tools/lor-transfer.cpp index 97695a07e8..607768b445 100644 --- a/miniapps/tools/lor-transfer.cpp +++ b/miniapps/tools/lor-transfer.cpp @@ -95,13 +95,7 @@ int main(int argc, char *argv[]) args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer", "-no-t", "--dont-use-pointwise-transfer", "Use pointwise transfer operators instead of L2 projection."); - args.Parse(); - if (!args.Good()) - { - args.PrintUsage(cout); - return 1; - } - args.PrintOptions(cout); + args.ParseCheck(); // Read the mesh from the given mesh file. Mesh mesh(mesh_file, 1, 1); @@ -157,6 +151,10 @@ int main(int argc, char *argv[]) direction = "HO -> LOR @ HO"; FunctionCoefficient RHO(RHO_exact); rho.ProjectCoefficient(RHO); + // Make sure AMR constraints are satisfied + rho.SetTrueVector(); + rho.SetFromTrueVector(); + double ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO "); if (vis) { visualize(HO_dc, "HO", Wx, Wy); Wx += offx; } @@ -193,17 +191,13 @@ int main(int argc, char *argv[]) } // HO* to LOR* dual fields - GridFunction ones(&fespace), ones_lor(&fespace_lor); - ones = 1.0; - ones_lor = 1.0; LinearForm M_rho(&fespace), M_rho_lor(&fespace_lor); if (!use_pointwise_transfer && gt->SupportsBackwardsOperator()) { const Operator &P = gt->BackwardOperator(); M_ho.Mult(rho, M_rho); P.MultTranspose(M_rho, M_rho_lor); - cout << "HO -> LOR dual field: " << fabs(M_rho(ones)-M_rho_lor(ones_lor)) - << endl << endl; + cout << "HO -> LOR dual field: " << abs(M_rho.Sum()-M_rho_lor.Sum()) << "\n\n"; } // LOR projections @@ -239,8 +233,7 @@ int main(int argc, char *argv[]) { M_lor.Mult(rho_lor, M_rho_lor); R.MultTranspose(M_rho_lor, M_rho); - cout << "LOR -> HO dual field: " << fabs(M_rho(ones)-M_rho_lor(ones_lor)) - << '\n'; + cout << "LOR -> HO dual field: " << abs(M_rho.Sum() - M_rho_lor.Sum()) << '\n'; } delete fec; @@ -288,14 +281,11 @@ double compute_mass(FiniteElementSpace *L2, double massL2, VisItDataCollection &dc, string prefix) { ConstantCoefficient one(1.0); - BilinearForm ML2(L2); - ML2.AddDomainIntegrator(new MassIntegrator(one)); - ML2.Assemble(); + LinearForm lf(L2); + lf.AddDomainIntegrator(new DomainLFIntegrator(one)); + lf.Assemble(); - GridFunction rhoone(L2); - rhoone = 1.0; - - double newmass = ML2.InnerProduct(*dc.GetField("density"),rhoone); + double newmass = lf(*dc.GetField("density")); cout.precision(18); cout << space << " " << prefix << " mass = " << newmass; if (massL2 >= 0) diff --git a/miniapps/tools/makefile b/miniapps/tools/makefile index e0ba96548a..196faceccb 100644 --- a/miniapps/tools/makefile +++ b/miniapps/tools/makefile @@ -28,7 +28,7 @@ MFEM_LIB_FILE = mfem_is_not_built SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer \ check-tmop-metric tmop-metric-magnitude -PAR_MINIAPPS = +PAR_MINIAPPS = plor-transfer ifeq ($(MFEM_USE_MPI),NO) MINIAPPS = $(SEQ_MINIAPPS) else @@ -81,9 +81,9 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP) @$(call mfem-test,$<,, Tools miniapp) # Testing: Specific execution options -# Do not test: display-basis, load-dc, convert-dc, get-values, lor-transfer +# Do not test: display-basis, load-dc, convert-dc, get-values, lor-transfer, plor-transfer NO_TEST_APPS = display-basis load-dc convert-dc get-values lor-transfer \ - check-tmop-metric tmop-metric-magnitude + plor-transfer check-tmop-metric tmop-metric-magnitude $(foreach app,$(NO_TEST_APPS),$(app)-test-seq $(app)-test-par): @true diff --git a/miniapps/tools/plor-transfer.cpp b/miniapps/tools/plor-transfer.cpp new file mode 100644 index 0000000000..d8d87b8e86 --- /dev/null +++ b/miniapps/tools/plor-transfer.cpp @@ -0,0 +1,360 @@ +// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced +// at the Lawrence Livermore National Laboratory. All Rights reserved. See files +// LICENSE and NOTICE for details. LLNL-CODE-806117. +// +// This file is part of the MFEM library. For more information and source code +// availability visit https://mfem.org. +// +// MFEM is free software; you can redistribute it and/or modify it under the +// terms of the BSD-3 license. We welcome feedback and contributions, see file +// CONTRIBUTING.md for details. +// +// ----------------------------------------------------------------------- +// Parallel LOR Transfer Miniapp: Map functions between HO and LOR spaces +// ----------------------------------------------------------------------- +// +// This miniapp visualizes the maps between a high-order (HO) finite element +// space, typically using high-order functions on a high-order mesh, and a +// low-order refined (LOR) finite element space, typically defined by 0th or 1st +// order functions on a low-order refinement of the HO mesh. +// +// The grid transfer operators are represented using either +// InterpolationGridTransfer or L2ProjectionGridTransfer (depending on the +// options requested by the user). The two transfer operators are then: +// +// 1. R: HO -> LOR, defined by GridTransfer::ForwardOperator +// 2. P: LOR -> HO, defined by GridTransfer::BackwardOperator +// +// While defined generally, these operators have some nice properties for +// particular finite element spaces. For example they satisfy PR=I, plus mass +// conservation in both directions for L2 fields. +// +// Compile with: make plor-transfer +// +// Sample runs: plor-transfer +// plor-transfer -h1 +// plor-transfer -t +// plor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4 +// plor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2 +// plor-transfer -lref 4 -o 4 -lo 0 -p 1 +// plor-transfer -lref 5 -o 4 -lo 0 -p 1 +// plor-transfer -lref 5 -o 4 -lo 3 -p 2 +// plor-transfer -lref 5 -o 4 -lo 0 -p 3 + +#include "mfem.hpp" +#include +#include + +using namespace std; +using namespace mfem; + +int problem = 1; // problem type + +int Wx = 0, Wy = 0; // window position +int Ww = 350, Wh = 350; // window size +int offx = Ww+5, offy = Wh+25; // window offsets + +string space; +string direction; + +// Exact functions to project +double RHO_exact(const Vector &x); + +// Helper functions +void visualize(VisItDataCollection &, string, int, int); +double compute_mass(ParFiniteElementSpace *, double, VisItDataCollection &, + string); + +int main(int argc, char *argv[]) +{ + // Initialize MPI and HYPRE. + Mpi::Init(argc, argv); + Hypre::Init(); + + // Parse command-line options. + const char *mesh_file = "../../data/star.mesh"; + int order = 3; + int lref = order+1; + int lorder = 0; + bool vis = true; + bool useH1 = false; + bool use_pointwise_transfer = false; + + OptionsParser args(argc, argv); + args.AddOption(&mesh_file, "-m", "--mesh", + "Mesh file to use."); + args.AddOption(&problem, "-p", "--problem", + "Problem type (see the RHO_exact function)."); + args.AddOption(&order, "-o", "--order", + "Finite element order (polynomial degree) or -1 for" + " isoparametric space."); + args.AddOption(&lref, "-lref", "--lor-ref-level", "LOR refinement level."); + args.AddOption(&lorder, "-lo", "--lor-order", + "LOR space order (polynomial degree, zero by default)."); + args.AddOption(&vis, "-vis", "--visualization", "-no-vis", + "--no-visualization", + "Enable or disable GLVis visualization."); + args.AddOption(&useH1, "-h1", "--use-h1", "-l2", "--use-l2", + "Use H1 spaces instead of L2."); + args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer", + "-no-t", "--dont-use-pointwise-transfer", + "Use pointwise transfer operators instead of L2 projection."); + args.ParseCheck(); + + // Read the mesh from the given mesh file. + Mesh serial_mesh(mesh_file, 1, 1); + ParMesh mesh(MPI_COMM_WORLD, serial_mesh); + serial_mesh.Clear(); + int dim = mesh.Dimension(); + + // Create the low-order refined mesh + int basis_lor = BasisType::GaussLobatto; // BasisType::ClosedUniform; + ParMesh mesh_lor = ParMesh::MakeRefined(mesh, lref, basis_lor); + + // Create spaces + FiniteElementCollection *fec, *fec_lor; + if (useH1) + { + space = "H1"; + if (lorder == 0) + { + lorder = 1; + if (Mpi::Root()) + { + cerr << "Switching the H1 LOR space order from 0 to 1\n"; + } + } + fec = new H1_FECollection(order, dim); + fec_lor = new H1_FECollection(lorder, dim); + } + else + { + space = "L2"; + fec = new L2_FECollection(order, dim); + fec_lor = new L2_FECollection(lorder, dim); + } + + ParFiniteElementSpace fespace(&mesh, fec); + ParFiniteElementSpace fespace_lor(&mesh_lor, fec_lor); + + ParGridFunction rho(&fespace); + ParGridFunction rho_lor(&fespace_lor); + + // Data collections for vis/analysis + VisItDataCollection HO_dc(MPI_COMM_WORLD, "HO", &mesh); + HO_dc.RegisterField("density", &rho); + VisItDataCollection LOR_dc(MPI_COMM_WORLD, "LOR", &mesh_lor); + LOR_dc.RegisterField("density", &rho_lor); + + ParBilinearForm M_ho(&fespace); + M_ho.AddDomainIntegrator(new MassIntegrator); + M_ho.Assemble(); + M_ho.Finalize(); + HypreParMatrix* M_ho_tdof = M_ho.ParallelAssemble(); + + ParBilinearForm M_lor(&fespace_lor); + M_lor.AddDomainIntegrator(new MassIntegrator); + M_lor.Assemble(); + M_lor.Finalize(); + HypreParMatrix* M_lor_tdof = M_lor.ParallelAssemble(); + + // HO projections + direction = "HO -> LOR @ HO"; + FunctionCoefficient RHO(RHO_exact); + rho.ProjectCoefficient(RHO); + // Make sure AMR constraints are satisfied + rho.SetTrueVector(); + rho.SetFromTrueVector(); + + double ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO "); + if (vis) { visualize(HO_dc, "HO", Wx, Wy); Wx += offx; } + + GridTransfer *gt; + if (use_pointwise_transfer) + { + gt = new InterpolationGridTransfer(fespace, fespace_lor); + } + else + { + gt = new L2ProjectionGridTransfer(fespace, fespace_lor); + } + const Operator &R = gt->ForwardOperator(); + + // HO->LOR restriction + direction = "HO -> LOR @ LOR"; + R.Mult(rho, rho_lor); + compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO) "); + if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy); Wx += offx; } + auto global_max = [](const Vector& v) + { + double max = v.Normlinf(); + MPI_Allreduce(MPI_IN_PLACE, &max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + return max; + }; + + if (gt->SupportsBackwardsOperator()) + { + const Operator &P = gt->BackwardOperator(); + // LOR->HO prolongation + direction = "HO -> LOR @ HO"; + ParGridFunction rho_prev = rho; + P.Mult(rho_lor, rho); + compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) "); + if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy); Wx = 0; Wy += offy; } + + rho_prev -= rho; + Vector rho_prev_true(fespace.GetTrueVSize()); + rho_prev.GetTrueDofs(rho_prev_true); + double l_inf = global_max(rho_prev_true); + if (Mpi::Root()) + { + cout.precision(12); + cout << "|HO - P(R(HO))|_∞ = " << l_inf << endl; + } + } + + // HO* to LOR* dual fields + ParLinearForm M_rho(&fespace), M_rho_lor(&fespace_lor); + auto global_sum = [](const Vector& v) + { + double sum = v.Sum(); + MPI_Allreduce(MPI_IN_PLACE, &sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + return sum; + }; + if (!use_pointwise_transfer && gt->SupportsBackwardsOperator()) + { + Vector M_rho_true(fespace.GetTrueVSize()); + M_ho_tdof->Mult(rho.GetTrueVector(), M_rho_true); + fespace.GetRestrictionOperator()->MultTranspose(M_rho_true, M_rho); + const Operator &P = gt->BackwardOperator(); + P.MultTranspose(M_rho, M_rho_lor); + double ho_dual_mass = global_sum(M_rho); + double lor_dual_mass = global_sum(M_rho_lor); + if (Mpi::Root()) + { + cout << "HO -> LOR dual field: " << abs(ho_dual_mass - lor_dual_mass) << "\n\n"; + } + } + + // LOR projections + direction = "LOR -> HO @ LOR"; + rho_lor.ProjectCoefficient(RHO); + ParGridFunction rho_lor_prev = rho_lor; + double lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR "); + if (vis) { visualize(LOR_dc, "LOR", Wx, Wy); Wx += offx; } + + if (gt->SupportsBackwardsOperator()) + { + const Operator &P = gt->BackwardOperator(); + // Prolongate to HO space + direction = "LOR -> HO @ HO"; + P.Mult(rho_lor, rho); + compute_mass(&fespace, lor_mass, HO_dc, "P(LOR) "); + if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy); Wx += offx; } + + // Restrict back to LOR space. This won't give the original function because + // the rho_lor doesn't necessarily live in the range of R. + direction = "LOR -> HO @ LOR"; + R.Mult(rho, rho_lor); + compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))"); + if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy); } + + rho_lor_prev -= rho_lor; + Vector rho_lor_prev_true(fespace_lor.GetTrueVSize()); + rho_lor_prev.GetTrueDofs(rho_lor_prev_true); + double l_inf = global_max(rho_lor_prev_true); + if (Mpi::Root()) + { + cout.precision(12); + cout << "|LOR - R(P(LOR))|_∞ = " << l_inf << endl; + } + } + + // LOR* to HO* dual fields + if (!use_pointwise_transfer) + { + Vector M_rho_lor_true(fespace_lor.GetTrueVSize()); + M_lor_tdof->Mult(rho_lor.GetTrueVector(), M_rho_lor_true); + fespace_lor.GetRestrictionOperator()->MultTranspose(M_rho_lor_true, + M_rho_lor); + R.MultTranspose(M_rho_lor, M_rho); + double ho_dual_mass = global_sum(M_rho); + double lor_dual_mass = global_sum(M_rho_lor); + + cout << lor_dual_mass << '\n'; + cout << ho_dual_mass << '\n'; + + if (Mpi::Root()) + { + cout << "LOR -> HO dual field: " << abs(ho_dual_mass - lor_dual_mass) << '\n'; + } + } + + delete fec; + delete fec_lor; + delete M_ho_tdof; + delete M_lor_tdof; + delete gt; + + return 0; +} + + +double RHO_exact(const Vector &x) +{ + switch (problem) + { + case 1: // smooth field + return x(1)+0.25*cos(2*M_PI*x.Norml2()); + case 2: // cubic function + return x(1)*x(1)*x(1) + 2*x(0)*x(1) + x(0); + case 3: // sharp gradient + return M_PI/2-atan(5*(2*x.Norml2()-1)); + case 4: // basis function + return (x.Norml2() < 0.1) ? 1 : 0; + default: + return 1.0; + } +} + + +void visualize(VisItDataCollection &dc, string prefix, int x, int y) +{ + int w = Ww, h = Wh; + + char vishost[] = "localhost"; + int visport = 19916; + + socketstream sol_sockL2(vishost, visport); + sol_sockL2 << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << + "\n"; + sol_sockL2.precision(8); + sol_sockL2 << "solution\n" << *dc.GetMesh() << *dc.GetField("density") + << "window_geometry " << x << " " << y << " " << w << " " << h + << "plot_caption '" << space << " " << prefix << " Density'" + << "window_title '" << direction << "'" << flush; +} + + +double compute_mass(ParFiniteElementSpace *L2, double massL2, + VisItDataCollection &dc, string prefix) +{ + ConstantCoefficient one(1.0); + ParLinearForm lf(L2); + lf.AddDomainIntegrator(new DomainLFIntegrator(one)); + lf.Assemble(); + + double newmass = lf(*dc.GetParField("density")); + if (Mpi::Root()) + { + cout.precision(18); + cout << space << " " << prefix << " mass = " << newmass; + if (massL2 >= 0) + { + cout.precision(4); + cout << " (" << fabs(newmass-massL2)*100/massL2 << "%)"; + } + cout << endl; + } + return newmass; +}