From b06b1c0006a84e33331dd153f51b990b776f3480 Mon Sep 17 00:00:00 2001 From: EB Chin Date: Mon, 20 Mar 2023 13:12:53 -0700 Subject: [PATCH 01/34] initial support for parallel h1 fields --- fem/transfer.cpp | 343 +++++++++++++++++++++++++++++------------------ fem/transfer.hpp | 47 +++++-- 2 files changed, 254 insertions(+), 136 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 40b34312fc..f5c91a13ac 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -12,6 +12,7 @@ #include "transfer.hpp" #include "bilinearform.hpp" #include "../general/forall.hpp" +#include "pfespace.hpp" namespace mfem { @@ -546,6 +547,121 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_) : L2Projection(fes_ho_, fes_lor_) { + // Basic PCG solver setup + pcg.SetPrintLevel(0); + pcg.SetMaxIter(1000); + // initial values for relative and absolute tolerance + SetRelTol(1e-13); + SetAbsTol(1e-13); +} + +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) + { + 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); + } +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( + 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_lor); + Vector y_dim(ndof_ho); + + for (int d = 0; d < vdim; ++d) + { + 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); + } +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( + 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_lor); + Vector y_dim(ndof_ho); + Vector xbar(ndof_ho); + + for (int d = 0; d < vdim; ++d) + { + 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); + } +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( + 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); + Vector xbar(ndof_ho); + + for (int d = 0; d < vdim; ++d) + { + 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); + } +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_) +{ + pcg.SetRelTol(p_rtol_); +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_) +{ + pcg.SetAbsTol(p_atol_); +} + +std::pair + L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() +{ + std::pair r_and_mlh; + Mesh* mesh_ho = fes_ho.GetMesh(); Mesh* mesh_lor = fes_lor.GetMesh(); int nel_ho = mesh_ho->GetNE(); @@ -553,7 +669,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(); @@ -617,12 +733,13 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( } // 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); + r_and_mlh.second = new SparseMatrix( + r_and_mlh.first->GetI(), r_and_mlh.first->GetJ(), NULL, + r_and_mlh.first->Height(), r_and_mlh.first->Width(), false, true, true); IntegrationPointTransformation ip_tr; IsoparametricTransformation& emb_tr = ip_tr.Transf; @@ -667,131 +784,15 @@ 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() -{ - 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) - { - 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); - } -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( - 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_lor); - Vector y_dim(ndof_ho); - - for (int d = 0; d < vdim; ++d) - { - 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); - } -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( - 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_lor); - Vector y_dim(ndof_ho); - Vector xbar(ndof_ho); - - for (int d = 0; d < vdim; ++d) - { - 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); - } -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( - 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); - Vector xbar(ndof_ho); - - for (int d = 0; d < vdim; ++d) - { - 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); - } -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_) -{ - pcg.SetRelTol(p_rtol_); -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_) -{ - pcg.SetAbsTol(p_atol_); -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() +SparseMatrix* L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() { const Table& elem_dof_ho = fes_ho.GetElementToDofTable(); const Table& elem_dof_lor = fes_lor.GetElementToDofTable(); @@ -871,13 +872,86 @@ 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; + SparseMatrix* 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::SerialL2ProjectionH1Space::SerialL2ProjectionH1Space( + const FiniteElementSpace& pfes_ho_, const FiniteElementSpace& pfes_lor_) +: L2ProjectionH1Space(pfes_ho_, pfes_lor_) +{ + std::tie(R_sm, M_LH_sm) = ComputeSparseRAndM_LH(); + RTxM_LH_sm = TransposeMult(*R_sm, *M_LH_sm); + + R = R_sm; + M_LH = M_LH_sm; + RTxM_LH = RTxM_LH_sm; + + Ds = DSmoother(*RTxM_LH_sm); + pcg.SetPreconditioner(Ds); + pcg.SetOperator(*RTxM_LH); +} + +L2ProjectionGridTransfer::SerialL2ProjectionH1Space::~SerialL2ProjectionH1Space() +{ + delete R_sm; + delete M_LH_sm; + delete RTxM_LH_sm; +} + +#ifdef MFEM_USE_MPI + +L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( + const ParFiniteElementSpace& pfes_ho_, const ParFiniteElementSpace& pfes_lor_) +: L2ProjectionH1Space(pfes_ho_, pfes_lor_) +{ + SparseMatrix* R_sm; + SparseMatrix* M_LH_sm; + std::tie(R_sm, M_LH_sm) = ComputeSparseRAndM_LH(); + SparseMatrix* RTxM_LH_sm = TransposeMult(*R_sm, *M_LH_sm); + + HypreParMatrix R_local = HypreParMatrix(pfes_ho_.GetComm(), + pfes_lor_.GlobalVSize(), pfes_ho_.GlobalVSize(), + pfes_lor_.GetDofOffsets(), pfes_ho_.GetDofOffsets(), R_sm); + HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho_.GetComm(), + pfes_lor_.GlobalVSize(), pfes_ho_.GlobalVSize(), + pfes_lor_.GetDofOffsets(), pfes_ho_.GetDofOffsets(), M_LH_sm); + HypreParMatrix RTxM_LH_local = HypreParMatrix(pfes_ho_.GetComm(), + pfes_ho_.GlobalVSize(), pfes_ho_.GetDofOffsets(), RTxM_LH_sm); + + delete R_sm; + delete M_LH_sm; + delete RTxM_LH_sm; + + R_par = mfem::RAP(pfes_lor_.Dof_TrueDof_Matrix(), &R_local, + pfes_ho_.Dof_TrueDof_Matrix()); + M_LH_par = mfem::RAP(pfes_lor_.Dof_TrueDof_Matrix(), &M_LH_local, + pfes_ho_.Dof_TrueDof_Matrix()); + RTxM_LH_par = mfem::RAP(&RTxM_LH_local, pfes_ho_.Dof_TrueDof_Matrix()); + + R = R_par; + M_LH = M_LH_par; + RTxM_LH = RTxM_LH_par; + + M = HypreBoomerAMG(*RTxM_LH_par); + pcg.SetPreconditioner(M); + pcg.SetOperator(*RTxM_LH); +} + +L2ProjectionGridTransfer::ParL2ProjectionH1Space::~ParL2ProjectionH1Space() +{ + delete R_par; + delete M_LH_par; + delete RTxM_LH_par; +} + +#endif + L2ProjectionGridTransfer::~L2ProjectionGridTransfer() { delete F; @@ -905,7 +979,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 SerialL2ProjectionH1Space(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 ParL2ProjectionH1Space(dom_pfes, ran_pfes); +#endif + } } else { diff --git a/fem/transfer.hpp b/fem/transfer.hpp index de937857ed..1911931465 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -258,20 +258,20 @@ protected: refined mesh (LOR). */ class L2ProjectionH1Space : public L2Projection { - // The restriction operator is represented as a SparseMatrix R. The + protected: + // 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. - SparseMatrix R; + Operator* R; // Used to compute P = (RTxM_LH)^(-1) M_LH^T - SparseMatrix M_LH; - SparseMatrix* RTxM_LH; + Operator* M_LH; + Operator* RTxM_LH; CGSolver pcg; - DSmoother Ds; - public: L2ProjectionH1Space(const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_); - virtual ~L2ProjectionH1Space(); + public: + virtual ~L2ProjectionH1Space() = default; /// 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,13 +305,44 @@ protected: virtual void ProlongateTranspose(const Vector& x, Vector& y) const; virtual void SetRelTol(double p_rtol_); virtual void SetAbsTol(double p_atol_); + protected: + std::pair ComputeSparseRAndM_LH(); private: /// Computes sparsity pattern and initializes R matrix. Based on /// BilinearForm::AllocMat() except maps between HO elements and LOR /// elements. - void AllocR(); + SparseMatrix* AllocR(); }; + class SerialL2ProjectionH1Space : public L2ProjectionH1Space + { + private: + SparseMatrix* R_sm; + SparseMatrix* M_LH_sm; + SparseMatrix* RTxM_LH_sm; + DSmoother Ds; + public: + SerialL2ProjectionH1Space(const FiniteElementSpace& pfes_ho_, + const FiniteElementSpace& pfes_lor_); + virtual ~SerialL2ProjectionH1Space(); + }; + + +#ifdef MFEM_USE_MPI + class ParL2ProjectionH1Space : public L2ProjectionH1Space + { + private: + HypreParMatrix* R_par; + HypreParMatrix* M_LH_par; + HypreParMatrix* RTxM_LH_par; + HypreBoomerAMG M; + public: + ParL2ProjectionH1Space(const ParFiniteElementSpace& pfes_ho_, + const ParFiniteElementSpace& pfes_lor_); + virtual ~ParL2ProjectionH1Space(); + }; +#endif + /** Mass-conservative prolongation operator going in the opposite direction as L2Projection. This operator is a left inverse to the L2Projection. */ class L2Prolongation : public Operator From 05f18737a0ab94c9a2a8205292335eb6263bbe73 Mon Sep 17 00:00:00 2001 From: EB Chin Date: Tue, 21 Mar 2023 09:18:43 -0700 Subject: [PATCH 02/34] initial example --- fem/transfer.cpp | 13 +- fem/transfer.hpp | 8 +- miniapps/tools/lor-transferp.cpp | 318 +++++++++++++++++++++++++++++++ 3 files changed, 330 insertions(+), 9 deletions(-) create mode 100644 miniapps/tools/lor-transferp.cpp diff --git a/fem/transfer.cpp b/fem/transfer.cpp index f5c91a13ac..4e2a709815 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -555,6 +555,13 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( SetAbsTol(1e-13); } +L2ProjectionGridTransfer::L2ProjectionH1Space::~L2ProjectionH1Space() +{ + delete R_sm; + delete M_LH_sm; + delete RTxM_LH_sm; +} + void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( const Vector& x, Vector& y) const { @@ -923,10 +930,6 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( pfes_lor_.GetDofOffsets(), pfes_ho_.GetDofOffsets(), M_LH_sm); HypreParMatrix RTxM_LH_local = HypreParMatrix(pfes_ho_.GetComm(), pfes_ho_.GlobalVSize(), pfes_ho_.GetDofOffsets(), RTxM_LH_sm); - - delete R_sm; - delete M_LH_sm; - delete RTxM_LH_sm; R_par = mfem::RAP(pfes_lor_.Dof_TrueDof_Matrix(), &R_local, pfes_ho_.Dof_TrueDof_Matrix()); @@ -939,8 +942,8 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( RTxM_LH = RTxM_LH_par; M = HypreBoomerAMG(*RTxM_LH_par); + pcg.SetOperator(*RTxM_LH_par); pcg.SetPreconditioner(M); - pcg.SetOperator(*RTxM_LH); } L2ProjectionGridTransfer::ParL2ProjectionH1Space::~ParL2ProjectionH1Space() diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 1911931465..98c855efd9 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -266,12 +266,15 @@ protected: // Used to compute P = (RTxM_LH)^(-1) M_LH^T Operator* M_LH; Operator* RTxM_LH; + SparseMatrix* R_sm; + SparseMatrix* M_LH_sm; + SparseMatrix* RTxM_LH_sm; CGSolver pcg; L2ProjectionH1Space(const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_); public: - virtual ~L2ProjectionH1Space() = default; + virtual ~L2ProjectionH1Space(); /// 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 @@ -317,9 +320,6 @@ protected: class SerialL2ProjectionH1Space : public L2ProjectionH1Space { private: - SparseMatrix* R_sm; - SparseMatrix* M_LH_sm; - SparseMatrix* RTxM_LH_sm; DSmoother Ds; public: SerialL2ProjectionH1Space(const FiniteElementSpace& pfes_ho_, diff --git a/miniapps/tools/lor-transferp.cpp b/miniapps/tools/lor-transferp.cpp new file mode 100644 index 0000000000..c8c75754ec --- /dev/null +++ b/miniapps/tools/lor-transferp.cpp @@ -0,0 +1,318 @@ +// Copyright (c) 2010-2022, 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. +// +// -------------------------------------------------------------- +// 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 lor-transfer +// +// Sample runs: lor-transfer +// lor-transfer -h1 +// lor-transfer -t +// lor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4 +// lor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2 +// lor-transfer -lref 4 -o 4 -lo 0 -p 1 +// lor-transfer -lref 5 -o 4 -lo 0 -p 1 +// lor-transfer -lref 5 -o 4 -lo 3 -p 2 +// lor-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.Parse(); + if (!args.Good()) + { + args.PrintUsage(cout); + return 1; + } + args.PrintOptions(cout); + + // 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; + 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(); + + ParBilinearForm M_lor(&fespace_lor); + M_lor.AddDomainIntegrator(new MassIntegrator); + M_lor.Assemble(); + M_lor.Finalize(); + + // HO projections + direction = "HO -> LOR @ HO"; + FunctionCoefficient RHO(RHO_exact); + rho.ProjectCoefficient(RHO); + 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; } + + 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; + cout.precision(12); + cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl; + } + + // HO* to LOR* dual fields + ParGridFunction ones(&fespace), ones_lor(&fespace_lor); + ones = 1.0; + ones_lor = 1.0; + ParLinearForm 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; + } + + // 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; + cout.precision(12); + cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl; + } + + // LOR* to HO* dual fields + if (!use_pointwise_transfer) + { + 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'; + } + + delete fec; + delete fec_lor; + + 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.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); + ParBilinearForm ML2(L2); + ML2.AddDomainIntegrator(new MassIntegrator(one)); + ML2.Assemble(); + ML2.Finalize(); + HypreParMatrix* pML2 = ML2.ParallelAssemble(); + + ParGridFunction rhoone(L2); + rhoone = 1.0; + + ParGridFunction Mdiag(L2); + pML2->Mult(rhoone, Mdiag); + double newmass = mfem::InnerProduct(*dc.GetParField("density"), Mdiag); + cout.precision(18); + cout << space << " " << prefix << " mass = " << newmass; + if (massL2 >= 0) + { + cout.precision(4); + cout << " (" << fabs(newmass-massL2)*100/massL2 << "%)"; + } + cout << endl; + return newmass; +} From a89d24843f69506af7dcaa003dc989641a0ef116 Mon Sep 17 00:00:00 2001 From: EB Chin Date: Thu, 23 Mar 2023 22:53:18 -0700 Subject: [PATCH 03/34] separate methods for parallel; example refinement --- fem/transfer.cpp | 402 +++++++++++++++++++++---------- fem/transfer.hpp | 82 +++++-- miniapps/tools/CMakeLists.txt | 5 + miniapps/tools/lor-transferp.cpp | 87 +++++-- 4 files changed, 407 insertions(+), 169 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 4e2a709815..c8d68590f5 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -546,122 +546,19 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::ProlongateTranspose( L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_) : L2Projection(fes_ho_, fes_lor_) -{ - // Basic PCG solver setup - pcg.SetPrintLevel(0); - pcg.SetMaxIter(1000); - // initial values for relative and absolute tolerance - SetRelTol(1e-13); - SetAbsTol(1e-13); -} +{} L2ProjectionGridTransfer::L2ProjectionH1Space::~L2ProjectionH1Space() -{ - delete R_sm; - delete M_LH_sm; - delete RTxM_LH_sm; -} - -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) - { - 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); - } -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( - 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_lor); - Vector y_dim(ndof_ho); - - for (int d = 0; d < vdim; ++d) - { - 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); - } -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( - 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_lor); - Vector y_dim(ndof_ho); - Vector xbar(ndof_ho); - - for (int d = 0; d < vdim; ++d) - { - 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); - } -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( - 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); - Vector xbar(ndof_ho); - - for (int d = 0; d < vdim; ++d) - { - 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); - } -} +{} void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_) { - pcg.SetRelTol(p_rtol_); + pcg.SetRelTol(sqrt(p_rtol_)); } void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_) { - pcg.SetAbsTol(p_atol_); + pcg.SetAbsTol(sqrt(p_atol_)); } std::pair @@ -734,10 +631,7 @@ std::pair } } // 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 r_and_mlh.first = AllocR(); @@ -899,6 +793,12 @@ L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SerialL2ProjectionH1Space( M_LH = M_LH_sm; RTxM_LH = RTxM_LH_sm; + // Basic PCG solver setup + pcg.SetPrintLevel(0); + pcg.SetMaxIter(1000); + // initial values for relative and absolute tolerance + SetRelTol(1e-20); + SetAbsTol(1e-20); Ds = DSmoother(*RTxM_LH_sm); pcg.SetPreconditioner(Ds); pcg.SetOperator(*RTxM_LH); @@ -911,39 +811,151 @@ L2ProjectionGridTransfer::SerialL2ProjectionH1Space::~SerialL2ProjectionH1Space( delete RTxM_LH_sm; } +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::LumpedMassInverse(Vector& ML_inv) const +{ + for (int i = 0; i < ML_inv.Size(); ++i) + { + ML_inv[i] = 1.0 / ML_inv[i]; + } +} + +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::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) + { + 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); + } +} + +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::MultTranspose( + 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_lor); + Vector y_dim(ndof_ho); + + for (int d = 0; d < vdim; ++d) + { + 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); + } +} + +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::Prolongate( + 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_lor); + Vector y_dim(ndof_ho); + Vector xbar(ndof_ho); + + for (int d = 0; d < vdim; ++d) + { + 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); + } +} + +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::ProlongateTranspose( + 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); + Vector xbar(ndof_ho); + + for (int d = 0; d < vdim; ++d) + { + 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); + } +} + #ifdef MFEM_USE_MPI L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( const ParFiniteElementSpace& pfes_ho_, const ParFiniteElementSpace& pfes_lor_) -: L2ProjectionH1Space(pfes_ho_, pfes_lor_) +: L2ProjectionH1Space(pfes_ho_, pfes_lor_), + pfes_ho(pfes_ho_), + pfes_lor(pfes_lor_), + pfes_ho_scalar(pfes_ho.GetParMesh(), pfes_ho.FEColl(), 1), + pfes_lor_scalar(pfes_lor.GetParMesh(), pfes_lor.FEColl(), 1) { SparseMatrix* R_sm; SparseMatrix* M_LH_sm; std::tie(R_sm, M_LH_sm) = ComputeSparseRAndM_LH(); - SparseMatrix* RTxM_LH_sm = TransposeMult(*R_sm, *M_LH_sm); - HypreParMatrix R_local = HypreParMatrix(pfes_ho_.GetComm(), - pfes_lor_.GlobalVSize(), pfes_ho_.GlobalVSize(), - pfes_lor_.GetDofOffsets(), pfes_ho_.GetDofOffsets(), R_sm); - HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho_.GetComm(), - pfes_lor_.GlobalVSize(), pfes_ho_.GlobalVSize(), - pfes_lor_.GetDofOffsets(), pfes_ho_.GetDofOffsets(), M_LH_sm); - HypreParMatrix RTxM_LH_local = HypreParMatrix(pfes_ho_.GetComm(), - pfes_ho_.GlobalVSize(), pfes_ho_.GetDofOffsets(), RTxM_LH_sm); + HypreParMatrix R_local = HypreParMatrix(pfes_ho_scalar.GetComm(), + pfes_lor_scalar.GlobalVSize(), pfes_ho_scalar.GlobalVSize(), + pfes_lor_scalar.GetDofOffsets(), pfes_ho_scalar.GetDofOffsets(), R_sm); + HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho_scalar.GetComm(), + pfes_lor_scalar.GlobalVSize(), pfes_ho_scalar.GlobalVSize(), + pfes_lor_scalar.GetDofOffsets(), pfes_ho_scalar.GetDofOffsets(), M_LH_sm); - R_par = mfem::RAP(pfes_lor_.Dof_TrueDof_Matrix(), &R_local, - pfes_ho_.Dof_TrueDof_Matrix()); - M_LH_par = mfem::RAP(pfes_lor_.Dof_TrueDof_Matrix(), &M_LH_local, - pfes_ho_.Dof_TrueDof_Matrix()); - RTxM_LH_par = mfem::RAP(&RTxM_LH_local, pfes_ho_.Dof_TrueDof_Matrix()); + R_par = RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), &R_local, + pfes_ho_scalar.Dof_TrueDof_Matrix()); + M_LH_par = RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), &M_LH_local, + pfes_ho_scalar.Dof_TrueDof_Matrix()); + HypreParMatrix* R_T = R_par->Transpose(); + RTxM_LH_par = ParMult(R_T, M_LH_par, true); R = R_par; M_LH = M_LH_par; RTxM_LH = RTxM_LH_par; + // Basic PCG solver setup + pcg = CGSolver(MPI_COMM_WORLD); + pcg.SetPrintLevel(0); + pcg.SetMaxIter(1000); + // initial values for relative and absolute tolerance + SetRelTol(1e-20); + SetAbsTol(1e-20); M = HypreBoomerAMG(*RTxM_LH_par); pcg.SetOperator(*RTxM_LH_par); - pcg.SetPreconditioner(M); + // pcg.SetPreconditioner(M); + + delete R_sm; + delete M_LH_sm; } L2ProjectionGridTransfer::ParL2ProjectionH1Space::~ParL2ProjectionH1Space() @@ -953,6 +965,142 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::~ParL2ProjectionH1Space() delete RTxM_LH_par; } +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::LumpedMassInverse(Vector& ML_inv) const +{ + Vector ML_inv_true(pfes_lor_scalar.GetTrueVSize()); + const Operator& P = *pfes_lor_scalar.GetProlongationMatrix(); + P.MultTranspose(ML_inv, ML_inv_true); + for (int i = 0; i < ML_inv_true.Size(); ++i) + { + ML_inv_true[i] = 1.0 / ML_inv_true[i]; + } + P.Mult(ML_inv_true, ML_inv); +} + +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::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); + const int ntdof_ho = pfes_ho_scalar.GetTrueVSize(); + const int ntdof_lor = pfes_lor_scalar.GetTrueVSize(); + Vector X_dim(ntdof_ho); + Vector Y_dim(ntdof_lor); + const Operator* fes_R = pfes_ho_scalar.GetRestrictionOperator(); + const Operator* fes_P = pfes_lor_scalar.GetProlongationMatrix(); + + for (int d = 0; d < vdim; ++d) + { + fes_ho.GetVDofs(d, dofs_ho); + fes_lor.GetVDofs(d, dofs_lor); + x.GetSubVector(dofs_ho, x_dim); + fes_R->Mult(x_dim, X_dim); + R->Mult(X_dim, Y_dim); + fes_P->Mult(Y_dim, y_dim); + y.SetSubVector(dofs_lor, y_dim); + } +} + +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::MultTranspose( + 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_lor); + Vector y_dim(ndof_ho); + const int ntdof_ho = pfes_ho_scalar.GetTrueVSize(); + const int ntdof_lor = pfes_lor_scalar.GetTrueVSize(); + Vector X_dim(ntdof_lor); + Vector Y_dim(ntdof_ho); + const Operator* fes_R = pfes_lor_scalar.GetRestrictionOperator(); + const Operator* fes_P = pfes_ho_scalar.GetProlongationMatrix(); + + for (int d = 0; d < vdim; ++d) + { + fes_ho.GetVDofs(d, dofs_ho); + fes_lor.GetVDofs(d, dofs_lor); + x.GetSubVector(dofs_lor, x_dim); + fes_R->Mult(x_dim, X_dim); + R->MultTranspose(X_dim, Y_dim); + fes_P->Mult(Y_dim, y_dim); + y.SetSubVector(dofs_ho, y_dim); + } +} + +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::Prolongate( + 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_lor); + Vector y_dim(ndof_ho); + const int ntdof_ho = pfes_ho_scalar.GetTrueVSize(); + const int ntdof_lor = pfes_lor_scalar.GetTrueVSize(); + Vector X_dim(ntdof_lor); + Vector Y_dim(ntdof_ho); + Vector Xbar(ntdof_ho); + const Operator* fes_R = pfes_lor_scalar.GetRestrictionOperator(); + const Operator* fes_P = pfes_ho_scalar.GetProlongationMatrix(); + + for (int d = 0; d < vdim; ++d) + { + fes_lor.GetVDofs(d, dofs_lor); + x.GetSubVector(dofs_lor, x_dim); + fes_R->Mult(x_dim, 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_P->Mult(Y_dim, y_dim); + fes_ho.GetVDofs(d, dofs_ho); + y.SetSubVector(dofs_ho, y_dim); + } +} + +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::ProlongateTranspose( + 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); + const int ntdof_ho = pfes_ho_scalar.GetTrueVSize(); + const int ntdof_lor = pfes_lor_scalar.GetTrueVSize(); + Vector X_dim(ntdof_ho); + Vector Y_dim(ntdof_lor); + const Operator* fes_R = pfes_ho_scalar.GetRestrictionOperator(); + const Operator* fes_P = pfes_lor_scalar.GetProlongationMatrix(); + Vector Xbar(ntdof_ho); + + for (int d = 0; d < vdim; ++d) + { + fes_ho.GetVDofs(d, dofs_ho); + x.GetSubVector(dofs_ho, x_dim); + fes_R->Mult(x_dim, 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_P->Mult(Y_dim, y_dim); + fes_lor.GetVDofs(d, dofs_lor); + y.SetSubVector(dofs_lor, y_dim); + } +} + #endif L2ProjectionGridTransfer::~L2ProjectionGridTransfer() diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 98c855efd9..018f1331be 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -266,15 +266,36 @@ protected: // Used to compute P = (RTxM_LH)^(-1) M_LH^T Operator* M_LH; Operator* RTxM_LH; - SparseMatrix* R_sm; - SparseMatrix* M_LH_sm; - SparseMatrix* RTxM_LH_sm; CGSolver pcg; L2ProjectionH1Space(const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_); public: virtual ~L2ProjectionH1Space(); + virtual void SetRelTol(double p_rtol_); + virtual void SetAbsTol(double p_atol_); + protected: + std::pair ComputeSparseRAndM_LH(); + virtual void LumpedMassInverse(Vector& ML_inv) const = 0; + private: + /// Computes sparsity pattern and initializes R matrix. Based on + /// BilinearForm::AllocMat() except maps between HO elements and LOR + /// elements. + SparseMatrix* AllocR(); + }; + + class SerialL2ProjectionH1Space : public L2ProjectionH1Space + { + private: + SparseMatrix* R_sm; + SparseMatrix* M_LH_sm; + SparseMatrix* RTxM_LH_sm; + DSmoother Ds; + virtual void LumpedMassInverse(Vector& ML_inv) const; + public: + SerialL2ProjectionH1Space(const FiniteElementSpace& pfes_ho_, + const FiniteElementSpace& pfes_lor_); + virtual ~SerialL2ProjectionH1Space(); /// 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 @@ -306,25 +327,6 @@ 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_); - protected: - std::pair ComputeSparseRAndM_LH(); - private: - /// Computes sparsity pattern and initializes R matrix. Based on - /// BilinearForm::AllocMat() except maps between HO elements and LOR - /// elements. - SparseMatrix* AllocR(); - }; - - class SerialL2ProjectionH1Space : public L2ProjectionH1Space - { - private: - DSmoother Ds; - public: - SerialL2ProjectionH1Space(const FiniteElementSpace& pfes_ho_, - const FiniteElementSpace& pfes_lor_); - virtual ~SerialL2ProjectionH1Space(); }; @@ -332,14 +334,50 @@ protected: class ParL2ProjectionH1Space : public L2ProjectionH1Space { private: + const ParFiniteElementSpace& pfes_ho; + const ParFiniteElementSpace& pfes_lor; + ParFiniteElementSpace pfes_ho_scalar; + ParFiniteElementSpace pfes_lor_scalar; HypreParMatrix* R_par; HypreParMatrix* M_LH_par; HypreParMatrix* RTxM_LH_par; HypreBoomerAMG M; + virtual void LumpedMassInverse(Vector& ML_inv) const; public: ParL2ProjectionH1Space(const ParFiniteElementSpace& pfes_ho_, const ParFiniteElementSpace& pfes_lor_); virtual ~ParL2ProjectionH1Space(); + /// 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 + /// finite element space. Refined mesh should be a uniform refinement of + /// the coarse mesh. Coefficients are computed through minimization of L2 + /// error between the fields. + virtual void Mult(const Vector& x, Vector& y) const; + /// Maps x, dual field coefficients defined on a refined mesh + /// with a low order H1 finite element space, to y, dual field + /// coefficients defined on a coarse mesh with a higher order H1 finite + /// element space. Refined mesh should be a uniform refinement of the + /// coarse mesh. Coefficients are computed through minimization of L2 + /// error between the primal fields. Note, if the x-coefficients + /// come from ProlongateTranspose, then mass is conserved. + virtual void MultTranspose(const Vector& x, Vector& y) const; + /// Maps x, primal field coefficients defined on a refined mesh + /// with a low order H1 finite element space, to y, primal field + /// coefficients defined on a coarse mesh with a higher order H1 finite + /// element space. Refined mesh should be a uniform refinement of the + /// coarse mesh. Coefficients are computed from the mass conservative + /// left-inverse prolongation operation. This functionality is also + /// provided as an Operator by L2Prolongation. + virtual void Prolongate(const Vector& x, Vector& y) const; + /// Maps x, dual field coefficients defined on a coarse mesh with + /// a higher order H1 finite element space, to y, dual field + /// coefficients defined on a refined mesh with a low order H1 finite + /// element space. Refined mesh should be a uniform refinement of the + /// coarse mesh. Coefficients are computed from the transpose of the mass + /// conservative left-inverse prolongation operation. This functionality + /// is also provided as an Operator by L2Prolongation. + virtual void ProlongateTranspose(const Vector& x, Vector& y) const; }; #endif diff --git a/miniapps/tools/CMakeLists.txt b/miniapps/tools/CMakeLists.txt index 7793e85a81..c2fd4efb07 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(lor-transferp + MAIN lor-transferp.cpp LIBRARIES mfem) +endif() + add_mfem_miniapp(check-tmop-metric MAIN check-tmop-metric.cpp LIBRARIES mfem) diff --git a/miniapps/tools/lor-transferp.cpp b/miniapps/tools/lor-transferp.cpp index c8c75754ec..dceaff4d1a 100644 --- a/miniapps/tools/lor-transferp.cpp +++ b/miniapps/tools/lor-transferp.cpp @@ -105,7 +105,11 @@ int main(int argc, char *argv[]) args.PrintUsage(cout); return 1; } - args.PrintOptions(cout); + if (Mpi::Root()) + { + args.PrintOptions(cout); + } + // Read the mesh from the given mesh file. Mesh serial_mesh(mesh_file, 1, 1); @@ -153,11 +157,13 @@ int main(int argc, char *argv[]) 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"; @@ -194,8 +200,11 @@ int main(int argc, char *argv[]) if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy); Wx = 0; Wy += offy; } rho_prev -= rho; - cout.precision(12); - cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl; + if (Mpi::Root()) + { + cout.precision(12); + cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl; + } } // HO* to LOR* dual fields @@ -206,10 +215,24 @@ int main(int argc, char *argv[]) if (!use_pointwise_transfer && gt->SupportsBackwardsOperator()) { const Operator &P = gt->BackwardOperator(); - M_ho.Mult(rho, M_rho); + Vector rho_true(rho.ParFESpace()->GetTrueVSize()); + rho.GetTrueDofs(rho_true); + Vector M_rho_true(M_rho.ParFESpace()->GetTrueVSize()); + M_ho_tdof->Mult(rho_true, M_rho_true); + M_rho.ParFESpace()->GetProlongationMatrix()->Mult(M_rho_true, 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; + Vector M_rho_lor_true(M_rho_lor.ParFESpace()->GetTrueVSize()); + M_rho_lor.ParFESpace()->GetRestrictionOperator()->Mult(M_rho_lor, M_rho_lor_true); + double local_ho_mass = M_rho_true.Sum(); + double local_lor_mass = M_rho_lor_true.Sum(); + double ho_mass; + double lor_mass; + MPI_Allreduce(&local_ho_mass, &ho_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&local_lor_mass, &lor_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + if (Mpi::Root()) + { + cout << "HO -> LOR dual field: " << fabs(ho_mass - lor_mass) << endl << endl; + } } // LOR projections @@ -236,21 +259,40 @@ int main(int argc, char *argv[]) if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy); } rho_lor_prev -= rho_lor; - cout.precision(12); - cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl; + if (Mpi::Root()) + { + cout.precision(12); + cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl; + } } // LOR* to HO* dual fields if (!use_pointwise_transfer) { - M_lor.Mult(rho_lor, M_rho_lor); + Vector rho_lor_true(rho_lor.ParFESpace()->GetTrueVSize()); + rho_lor.GetTrueDofs(rho_lor_true); + Vector M_rho_lor_true(M_rho_lor.ParFESpace()->GetTrueVSize()); + M_lor_tdof->Mult(rho_lor_true, M_rho_lor_true); + M_rho_lor.ParFESpace()->GetProlongationMatrix()->Mult(M_rho_lor_true, 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'; + Vector M_rho_true(M_rho.ParFESpace()->GetTrueVSize()); + M_rho.ParFESpace()->GetRestrictionOperator()->Mult(M_rho, M_rho_true); + double local_ho_mass = M_rho_true.Sum(); + double local_lor_mass = M_rho_lor_true.Sum(); + double ho_mass; + double lor_mass; + MPI_Allreduce(&local_ho_mass, &ho_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&local_lor_mass, &lor_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + if (Mpi::Root()) + { + cout << "LOR -> HO dual field: " << fabs(ho_mass - lor_mass) << '\n'; + } } delete fec; delete fec_lor; + delete M_ho_tdof; + delete M_lor_tdof; return 0; } @@ -300,19 +342,24 @@ double compute_mass(ParFiniteElementSpace *L2, double massL2, ML2.Finalize(); HypreParMatrix* pML2 = ML2.ParallelAssemble(); - ParGridFunction rhoone(L2); + Vector rhoone(L2->GetTrueVSize()); rhoone = 1.0; - ParGridFunction Mdiag(L2); + Vector Mdiag(L2->GetTrueVSize()); pML2->Mult(rhoone, Mdiag); - double newmass = mfem::InnerProduct(*dc.GetParField("density"), Mdiag); - cout.precision(18); - cout << space << " " << prefix << " mass = " << newmass; - if (massL2 >= 0) + HypreParVector* rho = dc.GetParField("density")->GetTrueDofs(); + double newmass = InnerProduct(MPI_COMM_WORLD, *rho, Mdiag); + delete rho; + if (Mpi::Root()) { - cout.precision(4); - cout << " (" << fabs(newmass-massL2)*100/massL2 << "%)"; + cout.precision(18); + cout << space << " " << prefix << " mass = " << newmass; + if (massL2 >= 0) + { + cout.precision(4); + cout << " (" << fabs(newmass-massL2)*100/massL2 << "%)"; + } + cout << endl; } - cout << endl; return newmass; } From ebdfac0dd78cd904e42b59a90a7d0d167c8f345a Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 24 Mar 2023 08:39:52 -0700 Subject: [PATCH 04/34] fix dual example --- miniapps/tools/lor-transferp.cpp | 47 ++++++++++++++++++++------------ 1 file changed, 30 insertions(+), 17 deletions(-) diff --git a/miniapps/tools/lor-transferp.cpp b/miniapps/tools/lor-transferp.cpp index dceaff4d1a..e6f0d33d55 100644 --- a/miniapps/tools/lor-transferp.cpp +++ b/miniapps/tools/lor-transferp.cpp @@ -9,9 +9,9 @@ // terms of the BSD-3 license. We welcome feedback and contributions, see file // CONTRIBUTING.md for details. // -// -------------------------------------------------------------- -// LOR Transfer Miniapp: Map functions between HO and LOR spaces -// -------------------------------------------------------------- +// ----------------------------------------------------------------------- +// 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 @@ -29,17 +29,17 @@ // particular finite element spaces. For example they satisfy PR=I, plus mass // conservation in both directions for L2 fields. // -// Compile with: make lor-transfer +// Compile with: make lor-transferp // -// Sample runs: lor-transfer -// lor-transfer -h1 -// lor-transfer -t -// lor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4 -// lor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2 -// lor-transfer -lref 4 -o 4 -lo 0 -p 1 -// lor-transfer -lref 5 -o 4 -lo 0 -p 1 -// lor-transfer -lref 5 -o 4 -lo 3 -p 2 -// lor-transfer -lref 5 -o 4 -lo 0 -p 3 +// Sample runs: lor-transferp +// lor-transferp -h1 +// lor-transferp -t +// lor-transferp -m ../../data/star-q2.mesh -lref 5 -p 4 +// lor-transferp -m ../../data/star-mixed.mesh -lref 3 -p 2 +// lor-transferp -lref 4 -o 4 -lo 0 -p 1 +// lor-transferp -lref 5 -o 4 -lo 0 -p 1 +// lor-transferp -lref 5 -o 4 -lo 3 -p 2 +// lor-transferp -lref 5 -o 4 -lo 0 -p 3 #include "mfem.hpp" #include @@ -102,14 +102,16 @@ int main(int argc, char *argv[]) args.Parse(); if (!args.Good()) { - args.PrintUsage(cout); + if (Mpi::Root()) + { + args.PrintUsage(cout); + } return 1; } if (Mpi::Root()) { args.PrintOptions(cout); } - // Read the mesh from the given mesh file. Mesh serial_mesh(mesh_file, 1, 1); @@ -200,10 +202,15 @@ int main(int argc, char *argv[]) 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_local = rho_prev_true.Normlinf(); + double l_inf; + MPI_Allreduce(&l_inf_local, &l_inf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); if (Mpi::Root()) { cout.precision(12); - cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl; + cout << "|HO - P(R(HO))|_∞ = " << l_inf << endl; } } @@ -259,10 +266,15 @@ int main(int argc, char *argv[]) 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_local = rho_lor_prev_true.Normlinf(); + double l_inf; + MPI_Allreduce(&l_inf_local, &l_inf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); if (Mpi::Root()) { cout.precision(12); - cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl; + cout << "|LOR - R(P(LOR))|_∞ = " << l_inf << endl; } } @@ -324,6 +336,7 @@ void visualize(VisItDataCollection &dc, string prefix, int x, int y) 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 From 16fede36224d89fb56b8b0936b05f0437a04d62e Mon Sep 17 00:00:00 2001 From: EB Chin Date: Fri, 24 Mar 2023 12:34:15 -0700 Subject: [PATCH 05/34] reorganize and consolidate parallel vs serial --- fem/transfer.cpp | 392 ++++++++++++++----------------- fem/transfer.hpp | 145 ++++++------ miniapps/tools/lor-transferp.cpp | 4 + 3 files changed, 261 insertions(+), 280 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index c8d68590f5..db917cf10a 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -551,6 +551,66 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( L2ProjectionGridTransfer::L2ProjectionH1Space::~L2ProjectionH1Space() {} +void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( + const Vector& x, Vector& y) const +{ + Vector X_dim(R->Width()); + Vector Y_dim(R->Height()); + for (int d = 0; d < fes_ho.GetVDim(); ++d) + { + GetHOTDofsByVDim(d, x, X_dim); + R->Mult(X_dim, Y_dim); + SetLORFromTDofsByVDim(d, Y_dim, y); + } +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( + const Vector& x, Vector& y) const +{ + Vector X_dim(R->Height()); + Vector Y_dim(R->Width()); + for (int d = 0; d < fes_ho.GetVDim(); ++d) + { + GetLORTDofsByVDim(d, x, X_dim); + R->MultTranspose(X_dim, Y_dim); + SetHOFromTDofsByVDim(d, Y_dim, y); + } +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( + const Vector& x, Vector& y) const +{ + Vector X_dim(M_LH->Height()); + Vector Y_dim(pcg.Height()); + Vector Xbar(pcg.Width()); + for (int d = 0; d < fes_ho.GetVDim(); ++d) + { + GetLORTDofsByVDim(d, x, 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); + SetHOFromTDofsByVDim(d, Y_dim, y); + } +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( + const Vector& x, Vector& y) const +{ + Vector X_dim(pcg.Width()); + Vector Y_dim(M_LH->Height()); + Vector Xbar(pcg.Height()); + for (int d = 0; d < fes_ho.GetVDim(); ++d) + { + GetHOTDofsByVDim(d, x, 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); + SetLORFromTDofsByVDim(d, Y_dim, y); + } +} + void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_) { pcg.SetRelTol(sqrt(p_rtol_)); @@ -783,8 +843,8 @@ SparseMatrix* L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() } L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SerialL2ProjectionH1Space( - const FiniteElementSpace& pfes_ho_, const FiniteElementSpace& pfes_lor_) -: L2ProjectionH1Space(pfes_ho_, pfes_lor_) + const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_) +: L2ProjectionH1Space(fes_ho_, fes_lor_) { std::tie(R_sm, M_LH_sm) = ComputeSparseRAndM_LH(); RTxM_LH_sm = TransposeMult(*R_sm, *M_LH_sm); @@ -819,96 +879,54 @@ void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::LumpedMassInverse(Vect } } -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::Mult( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::GetHOTDofsByVDim( + int vdim, const Vector& x, Vector& x_vdim_true) 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) - { - 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); - } + GetDofsByVDim(fes_ho, vdim, x, x_vdim_true); } -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::MultTranspose( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::GetLORTDofsByVDim( + int vdim, const Vector& x, Vector& x_vdim_true) 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) - { - 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); - } + GetDofsByVDim(fes_lor, vdim, x, x_vdim_true); } -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::Prolongate( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::GetDofsByVDim( + const FiniteElementSpace& fes, + int vdim, const Vector& x, Vector& x_vdim) 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) - { - 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); - } + MFEM_ASSERT(x.Size() == fes.GetVSize(), + "Size of x Vector must match the number of vector DOFs."); + MFEM_ASSERT(x_vdim.Size() == fes.GetNDofs(), + "Size of x_vdim Vector must match the number of DOFs."); + Array vdim_dofs(fes.GetNDofs()); + fes.GetVDofs(vdim, vdim_dofs); + x.GetSubVector(vdim_dofs, x_vdim); } -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::ProlongateTranspose( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SetHOFromTDofsByVDim( + int vdim, const Vector& y_vdim_true, 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); - Vector xbar(ndof_ho); + SetFromDofsByVDim(fes_ho, vdim, y_vdim_true, y); +} - for (int d = 0; d < vdim; ++d) - { - 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); - } +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SetLORFromTDofsByVDim( + int vdim, const Vector& y_vdim_true, Vector& y) const +{ + SetFromDofsByVDim(fes_lor, vdim, y_vdim_true, y); +} + +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SetFromDofsByVDim( + const FiniteElementSpace& fes, + int vdim, const Vector& y_vdim, Vector& y) const +{ + MFEM_ASSERT(y_vdim.Size() == fes.GetNDofs(), + "Size of y_vdim Vector must match the number of DOFs."); + MFEM_ASSERT(y.Size() == fes.GetVSize(), + "Size of y Vector must match the number of vector DOFs."); + Array vdim_dofs(fes.GetNDofs()); + fes.GetVDofs(vdim, vdim_dofs); + y.SetSubVector(vdim_dofs, y_vdim); } #ifdef MFEM_USE_MPI @@ -917,25 +935,35 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( const ParFiniteElementSpace& pfes_ho_, const ParFiniteElementSpace& pfes_lor_) : L2ProjectionH1Space(pfes_ho_, pfes_lor_), pfes_ho(pfes_ho_), - pfes_lor(pfes_lor_), - pfes_ho_scalar(pfes_ho.GetParMesh(), pfes_ho.FEColl(), 1), - pfes_lor_scalar(pfes_lor.GetParMesh(), pfes_lor.FEColl(), 1) + pfes_lor(pfes_lor_) { + if (pfes_ho.GetVDim() != 1) + { + pfes_ho_scalar = new ParFiniteElementSpace(pfes_ho.GetParMesh(), + pfes_ho.FEColl(), 1); + pfes_lor_scalar = new ParFiniteElementSpace(pfes_lor.GetParMesh(), + pfes_lor.FEColl(), 1); + } + else // vdim == 1 + { + pfes_ho_scalar = &pfes_ho; + pfes_lor_scalar = &pfes_lor; + } SparseMatrix* R_sm; SparseMatrix* M_LH_sm; std::tie(R_sm, M_LH_sm) = ComputeSparseRAndM_LH(); - HypreParMatrix R_local = HypreParMatrix(pfes_ho_scalar.GetComm(), - pfes_lor_scalar.GlobalVSize(), pfes_ho_scalar.GlobalVSize(), - pfes_lor_scalar.GetDofOffsets(), pfes_ho_scalar.GetDofOffsets(), R_sm); - HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho_scalar.GetComm(), - pfes_lor_scalar.GlobalVSize(), pfes_ho_scalar.GlobalVSize(), - pfes_lor_scalar.GetDofOffsets(), pfes_ho_scalar.GetDofOffsets(), M_LH_sm); + HypreParMatrix R_local = HypreParMatrix(pfes_ho_scalar->GetComm(), + pfes_lor_scalar->GlobalVSize(), pfes_ho_scalar->GlobalVSize(), + pfes_lor_scalar->GetDofOffsets(), pfes_ho_scalar->GetDofOffsets(), R_sm); + HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho_scalar->GetComm(), + pfes_lor_scalar->GlobalVSize(), pfes_ho_scalar->GlobalVSize(), + pfes_lor_scalar->GetDofOffsets(), pfes_ho_scalar->GetDofOffsets(), M_LH_sm); - R_par = RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), &R_local, - pfes_ho_scalar.Dof_TrueDof_Matrix()); - M_LH_par = RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), &M_LH_local, - pfes_ho_scalar.Dof_TrueDof_Matrix()); + R_par = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(), &R_local, + pfes_ho_scalar->Dof_TrueDof_Matrix()); + M_LH_par = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(), &M_LH_local, + pfes_ho_scalar->Dof_TrueDof_Matrix()); HypreParMatrix* R_T = R_par->Transpose(); RTxM_LH_par = ParMult(R_T, M_LH_par, true); @@ -948,27 +976,35 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( pcg.SetPrintLevel(0); pcg.SetMaxIter(1000); // initial values for relative and absolute tolerance - SetRelTol(1e-20); - SetAbsTol(1e-20); - M = HypreBoomerAMG(*RTxM_LH_par); + SetRelTol(1e-25); + SetAbsTol(1e-25); + M = new HypreBoomerAMG(*RTxM_LH_par); + M->SetPrintLevel(0); + pcg.SetPreconditioner(*M); pcg.SetOperator(*RTxM_LH_par); - // pcg.SetPreconditioner(M); delete R_sm; delete M_LH_sm; + delete R_T; } L2ProjectionGridTransfer::ParL2ProjectionH1Space::~ParL2ProjectionH1Space() { + if (pfes_ho.GetVDim() != 1) + { + delete pfes_ho_scalar; + delete pfes_lor_scalar; + } delete R_par; delete M_LH_par; delete RTxM_LH_par; + delete M; } void L2ProjectionGridTransfer::ParL2ProjectionH1Space::LumpedMassInverse(Vector& ML_inv) const { - Vector ML_inv_true(pfes_lor_scalar.GetTrueVSize()); - const Operator& P = *pfes_lor_scalar.GetProlongationMatrix(); + Vector ML_inv_true(pfes_lor_scalar->GetTrueVSize()); + const Operator& P = *pfes_lor_scalar->GetProlongationMatrix(); P.MultTranspose(ML_inv, ML_inv_true); for (int i = 0; i < ML_inv_true.Size(); ++i) { @@ -977,128 +1013,62 @@ void L2ProjectionGridTransfer::ParL2ProjectionH1Space::LumpedMassInverse(Vector& P.Mult(ML_inv_true, ML_inv); } -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::Mult( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::GetHOTDofsByVDim( + int vdim, const Vector& x, Vector& x_vdim_true) 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); - const int ntdof_ho = pfes_ho_scalar.GetTrueVSize(); - const int ntdof_lor = pfes_lor_scalar.GetTrueVSize(); - Vector X_dim(ntdof_ho); - Vector Y_dim(ntdof_lor); - const Operator* fes_R = pfes_ho_scalar.GetRestrictionOperator(); - const Operator* fes_P = pfes_lor_scalar.GetProlongationMatrix(); - - for (int d = 0; d < vdim; ++d) - { - fes_ho.GetVDofs(d, dofs_ho); - fes_lor.GetVDofs(d, dofs_lor); - x.GetSubVector(dofs_ho, x_dim); - fes_R->Mult(x_dim, X_dim); - R->Mult(X_dim, Y_dim); - fes_P->Mult(Y_dim, y_dim); - y.SetSubVector(dofs_lor, y_dim); - } + GetTDofsByVDim(pfes_ho, *pfes_ho_scalar, vdim, x, x_vdim_true); } -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::MultTranspose( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::GetLORTDofsByVDim( + int vdim, const Vector& x, Vector& x_vdim_true) 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); - const int ntdof_ho = pfes_ho_scalar.GetTrueVSize(); - const int ntdof_lor = pfes_lor_scalar.GetTrueVSize(); - Vector X_dim(ntdof_lor); - Vector Y_dim(ntdof_ho); - const Operator* fes_R = pfes_lor_scalar.GetRestrictionOperator(); - const Operator* fes_P = pfes_ho_scalar.GetProlongationMatrix(); - - for (int d = 0; d < vdim; ++d) - { - fes_ho.GetVDofs(d, dofs_ho); - fes_lor.GetVDofs(d, dofs_lor); - x.GetSubVector(dofs_lor, x_dim); - fes_R->Mult(x_dim, X_dim); - R->MultTranspose(X_dim, Y_dim); - fes_P->Mult(Y_dim, y_dim); - y.SetSubVector(dofs_ho, y_dim); - } + GetTDofsByVDim(pfes_lor, *pfes_lor_scalar, vdim, x, x_vdim_true); } -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::Prolongate( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::GetTDofsByVDim( + const ParFiniteElementSpace& pfes, const ParFiniteElementSpace& pfes_scalar, + int vdim, const Vector& x, Vector& x_vdim_true) 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); - const int ntdof_ho = pfes_ho_scalar.GetTrueVSize(); - const int ntdof_lor = pfes_lor_scalar.GetTrueVSize(); - Vector X_dim(ntdof_lor); - Vector Y_dim(ntdof_ho); - Vector Xbar(ntdof_ho); - const Operator* fes_R = pfes_lor_scalar.GetRestrictionOperator(); - const Operator* fes_P = pfes_ho_scalar.GetProlongationMatrix(); - - for (int d = 0; d < vdim; ++d) - { - fes_lor.GetVDofs(d, dofs_lor); - x.GetSubVector(dofs_lor, x_dim); - fes_R->Mult(x_dim, 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_P->Mult(Y_dim, y_dim); - fes_ho.GetVDofs(d, dofs_ho); - y.SetSubVector(dofs_ho, y_dim); - } + MFEM_ASSERT(x.Size() == pfes.GetVSize(), + "Size of x Vector must match the number of vector DOFs."); + MFEM_ASSERT(x_vdim_true.Size() == pfes_scalar.GetTrueVSize(), + "Size of x_vdim_true Vector must match the number of scalar true DOFs."); + // transfer to vector of vdim dofs + Array vdim_dofs(pfes.GetNDofs()); + pfes.GetVDofs(vdim, vdim_dofs); + Vector x_vdim(pfes.GetNDofs()); + x.GetSubVector(vdim_dofs, x_vdim); + // transfer to vector of true dofs on the vdim + pfes_scalar.GetRestrictionOperator()->Mult(x_vdim, x_vdim_true); } -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::ProlongateTranspose( - const Vector& x, Vector& y) const +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::SetHOFromTDofsByVDim( + int vdim, const Vector& y_vdim_true, 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); - const int ntdof_ho = pfes_ho_scalar.GetTrueVSize(); - const int ntdof_lor = pfes_lor_scalar.GetTrueVSize(); - Vector X_dim(ntdof_ho); - Vector Y_dim(ntdof_lor); - const Operator* fes_R = pfes_ho_scalar.GetRestrictionOperator(); - const Operator* fes_P = pfes_lor_scalar.GetProlongationMatrix(); - Vector Xbar(ntdof_ho); + SetFromTDofsByVDim(pfes_ho, *pfes_ho_scalar, vdim, y_vdim_true, y); +} - for (int d = 0; d < vdim; ++d) - { - fes_ho.GetVDofs(d, dofs_ho); - x.GetSubVector(dofs_ho, x_dim); - fes_R->Mult(x_dim, 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_P->Mult(Y_dim, y_dim); - fes_lor.GetVDofs(d, dofs_lor); - y.SetSubVector(dofs_lor, y_dim); - } +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::SetLORFromTDofsByVDim( + int vdim, const Vector& y_vdim_true, Vector& y) const +{ + SetFromTDofsByVDim(pfes_lor, *pfes_lor_scalar, vdim, y_vdim_true, y); +} + +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::SetFromTDofsByVDim( + const ParFiniteElementSpace& pfes, const ParFiniteElementSpace& pfes_scalar, + int vdim, const Vector& y_vdim_true, Vector& y) const +{ + MFEM_ASSERT(y_vdim_true.Size() == pfes_scalar.GetTrueVSize(), + "Size of y_vdim_true Vector must match the number of scalar true DOFs."); + MFEM_ASSERT(y.Size() == pfes.GetVSize(), + "Size of y Vector must match the number of vector DOFs."); + // prologate to vector of dofs on the vdim + Vector y_vdim(pfes.GetNDofs()); + pfes_scalar.GetProlongationMatrix()->Mult(y_vdim_true, y_vdim); + // transfer to global vector + Array vdim_dofs(pfes.GetNDofs()); + pfes.GetVDofs(vdim, vdim_dofs); + y.SetSubVector(vdim_dofs, y_vdim); } #endif diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 018f1331be..2cf0ffba1a 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -272,11 +272,50 @@ protected: const FiniteElementSpace& fes_lor_); public: virtual ~L2ProjectionH1Space(); + /// 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 + /// finite element space. Refined mesh should be a uniform refinement of + /// the coarse mesh. Coefficients are computed through minimization of L2 + /// error between the fields. + virtual void Mult(const Vector& x, Vector& y) const; + /// Maps x, dual field coefficients defined on a refined mesh + /// with a low order H1 finite element space, to y, dual field + /// coefficients defined on a coarse mesh with a higher order H1 finite + /// element space. Refined mesh should be a uniform refinement of the + /// coarse mesh. Coefficients are computed through minimization of L2 + /// error between the primal fields. Note, if the x-coefficients + /// come from ProlongateTranspose, then mass is conserved. + virtual void MultTranspose(const Vector& x, Vector& y) const; + /// Maps x, primal field coefficients defined on a refined mesh + /// with a low order H1 finite element space, to y, primal field + /// coefficients defined on a coarse mesh with a higher order H1 finite + /// element space. Refined mesh should be a uniform refinement of the + /// coarse mesh. Coefficients are computed from the mass conservative + /// left-inverse prolongation operation. This functionality is also + /// provided as an Operator by L2Prolongation. + virtual void Prolongate(const Vector& x, Vector& y) const; + /// Maps x, dual field coefficients defined on a coarse mesh with + /// a higher order H1 finite element space, to y, dual field + /// coefficients defined on a refined mesh with a low order H1 finite + /// element space. Refined mesh should be a uniform refinement of the + /// coarse mesh. Coefficients are computed from the transpose of the mass + /// 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_); protected: std::pair ComputeSparseRAndM_LH(); virtual void LumpedMassInverse(Vector& ML_inv) const = 0; + virtual void GetHOTDofsByVDim(int vdim, const Vector& x, + Vector& x_vdim_true) const = 0; + virtual void GetLORTDofsByVDim(int vdim, const Vector& x, + Vector& x_vdim_true) const = 0; + virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + Vector& y) const = 0; + virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + Vector& y) const = 0; private: /// Computes sparsity pattern and initializes R matrix. Based on /// BilinearForm::AllocMat() except maps between HO elements and LOR @@ -291,42 +330,25 @@ protected: SparseMatrix* M_LH_sm; SparseMatrix* RTxM_LH_sm; DSmoother Ds; - virtual void LumpedMassInverse(Vector& ML_inv) const; public: - SerialL2ProjectionH1Space(const FiniteElementSpace& pfes_ho_, - const FiniteElementSpace& pfes_lor_); + SerialL2ProjectionH1Space(const FiniteElementSpace& fes_ho_, + const FiniteElementSpace& fes_lor_); virtual ~SerialL2ProjectionH1Space(); - /// 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 - /// finite element space. Refined mesh should be a uniform refinement of - /// the coarse mesh. Coefficients are computed through minimization of L2 - /// error between the fields. - virtual void Mult(const Vector& x, Vector& y) const; - /// Maps x, dual field coefficients defined on a refined mesh - /// with a low order H1 finite element space, to y, dual field - /// coefficients defined on a coarse mesh with a higher order H1 finite - /// element space. Refined mesh should be a uniform refinement of the - /// coarse mesh. Coefficients are computed through minimization of L2 - /// error between the primal fields. Note, if the x-coefficients - /// come from ProlongateTranspose, then mass is conserved. - virtual void MultTranspose(const Vector& x, Vector& y) const; - /// Maps x, primal field coefficients defined on a refined mesh - /// with a low order H1 finite element space, to y, primal field - /// coefficients defined on a coarse mesh with a higher order H1 finite - /// element space. Refined mesh should be a uniform refinement of the - /// coarse mesh. Coefficients are computed from the mass conservative - /// left-inverse prolongation operation. This functionality is also - /// provided as an Operator by L2Prolongation. - virtual void Prolongate(const Vector& x, Vector& y) const; - /// Maps x, dual field coefficients defined on a coarse mesh with - /// a higher order H1 finite element space, to y, dual field - /// coefficients defined on a refined mesh with a low order H1 finite - /// element space. Refined mesh should be a uniform refinement of the - /// coarse mesh. Coefficients are computed from the transpose of the mass - /// conservative left-inverse prolongation operation. This functionality - /// is also provided as an Operator by L2Prolongation. - virtual void ProlongateTranspose(const Vector& x, Vector& y) const; + private: + /// Computes inverse of a lumped mass matrix (stored as a vector) + virtual void LumpedMassInverse(Vector& ML_inv) const; + virtual void GetHOTDofsByVDim(int vdim, const Vector& x, + Vector& x_vdim_true) const; + virtual void GetLORTDofsByVDim(int vdim, const Vector& x, + Vector& x_vdim_true) const; + virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + Vector& y) const; + virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + Vector& y) const; + void GetDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& x, Vector& x_vdim) const; + void SetFromDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& y_vdim, Vector& y) const; }; @@ -336,48 +358,33 @@ protected: private: const ParFiniteElementSpace& pfes_ho; const ParFiniteElementSpace& pfes_lor; - ParFiniteElementSpace pfes_ho_scalar; - ParFiniteElementSpace pfes_lor_scalar; + const ParFiniteElementSpace* pfes_ho_scalar; + const ParFiniteElementSpace* pfes_lor_scalar; HypreParMatrix* R_par; HypreParMatrix* M_LH_par; HypreParMatrix* RTxM_LH_par; - HypreBoomerAMG M; - virtual void LumpedMassInverse(Vector& ML_inv) const; + HypreBoomerAMG* M; public: ParL2ProjectionH1Space(const ParFiniteElementSpace& pfes_ho_, const ParFiniteElementSpace& pfes_lor_); virtual ~ParL2ProjectionH1Space(); - /// 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 - /// finite element space. Refined mesh should be a uniform refinement of - /// the coarse mesh. Coefficients are computed through minimization of L2 - /// error between the fields. - virtual void Mult(const Vector& x, Vector& y) const; - /// Maps x, dual field coefficients defined on a refined mesh - /// with a low order H1 finite element space, to y, dual field - /// coefficients defined on a coarse mesh with a higher order H1 finite - /// element space. Refined mesh should be a uniform refinement of the - /// coarse mesh. Coefficients are computed through minimization of L2 - /// error between the primal fields. Note, if the x-coefficients - /// come from ProlongateTranspose, then mass is conserved. - virtual void MultTranspose(const Vector& x, Vector& y) const; - /// Maps x, primal field coefficients defined on a refined mesh - /// with a low order H1 finite element space, to y, primal field - /// coefficients defined on a coarse mesh with a higher order H1 finite - /// element space. Refined mesh should be a uniform refinement of the - /// coarse mesh. Coefficients are computed from the mass conservative - /// left-inverse prolongation operation. This functionality is also - /// provided as an Operator by L2Prolongation. - virtual void Prolongate(const Vector& x, Vector& y) const; - /// Maps x, dual field coefficients defined on a coarse mesh with - /// a higher order H1 finite element space, to y, dual field - /// coefficients defined on a refined mesh with a low order H1 finite - /// element space. Refined mesh should be a uniform refinement of the - /// coarse mesh. Coefficients are computed from the transpose of the mass - /// conservative left-inverse prolongation operation. This functionality - /// is also provided as an Operator by L2Prolongation. - virtual void ProlongateTranspose(const Vector& x, Vector& y) const; + private: + /// Computes inverse of a lumped mass matrix (stored as a vector) + virtual void LumpedMassInverse(Vector& ML_inv) const; + virtual void GetHOTDofsByVDim(int vdim, const Vector& x, + Vector& x_vdim_true) const; + virtual void GetLORTDofsByVDim(int vdim, const Vector& x, + Vector& x_vdim_true) const; + virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + Vector& y) const; + virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + Vector& y) const; + void GetTDofsByVDim(const ParFiniteElementSpace& pfes, + const ParFiniteElementSpace& pfes_scalar, + int vdim, const Vector& x, Vector& x_vdim_true) const; + void SetFromTDofsByVDim(const ParFiniteElementSpace& pfes, + const ParFiniteElementSpace& pfes_scalar, + int vdim, const Vector& y_vdim_true, Vector& y) const; }; #endif diff --git a/miniapps/tools/lor-transferp.cpp b/miniapps/tools/lor-transferp.cpp index e6f0d33d55..abcd885cb3 100644 --- a/miniapps/tools/lor-transferp.cpp +++ b/miniapps/tools/lor-transferp.cpp @@ -301,10 +301,13 @@ int main(int argc, char *argv[]) } } + Mpi::Finalize(); + delete fec; delete fec_lor; delete M_ho_tdof; delete M_lor_tdof; + delete gt; return 0; } @@ -360,6 +363,7 @@ double compute_mass(ParFiniteElementSpace *L2, double massL2, Vector Mdiag(L2->GetTrueVSize()); pML2->Mult(rhoone, Mdiag); + delete pML2; HypreParVector* rho = dc.GetParField("density")->GetTrueDofs(); double newmass = InnerProduct(MPI_COMM_WORLD, *rho, Mdiag); delete rho; From ada5c8f4e9b17f3e67100217cada1b89a0776673 Mon Sep 17 00:00:00 2001 From: EB Chin Date: Fri, 24 Mar 2023 12:55:34 -0700 Subject: [PATCH 06/34] update docs --- fem/transfer.cpp | 24 +++++++++++------------ fem/transfer.hpp | 51 +++++++++++++++++++++++++++++++++++++++++++----- 2 files changed, 58 insertions(+), 17 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index db917cf10a..d0c1e52632 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -621,6 +621,16 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_) pcg.SetAbsTol(sqrt(p_atol_)); } +void L2ProjectionGridTransfer::L2ProjectionH1Space::InitializeCGSolver() +{ + // Basic PCG solver setup + pcg.SetPrintLevel(0); + pcg.SetMaxIter(1000); + // initial values for relative and absolute tolerance + SetRelTol(1e-25); + SetAbsTol(1e-25); +} + std::pair L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() { @@ -853,12 +863,7 @@ L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SerialL2ProjectionH1Space( M_LH = M_LH_sm; RTxM_LH = RTxM_LH_sm; - // Basic PCG solver setup - pcg.SetPrintLevel(0); - pcg.SetMaxIter(1000); - // initial values for relative and absolute tolerance - SetRelTol(1e-20); - SetAbsTol(1e-20); + InitializeCGSolver(); Ds = DSmoother(*RTxM_LH_sm); pcg.SetPreconditioner(Ds); pcg.SetOperator(*RTxM_LH); @@ -971,13 +976,8 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( M_LH = M_LH_par; RTxM_LH = RTxM_LH_par; - // Basic PCG solver setup pcg = CGSolver(MPI_COMM_WORLD); - pcg.SetPrintLevel(0); - pcg.SetMaxIter(1000); - // initial values for relative and absolute tolerance - SetRelTol(1e-25); - SetAbsTol(1e-25); + InitializeCGSolver(); M = new HypreBoomerAMG(*RTxM_LH_par); M->SetPrintLevel(0); pcg.SetPreconditioner(*M); diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 2cf0ffba1a..4f31c59c3b 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -253,9 +253,9 @@ protected: virtual void SetAbsTol(double p_atol_) {} }; - /** 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). */ + /** Abstract 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). */ class L2ProjectionH1Space : public L2Projection { protected: @@ -306,23 +306,40 @@ protected: virtual void SetRelTol(double p_rtol_); virtual void SetAbsTol(double p_atol_); protected: + /// Sets default shared options for the CGSolver + void InitializeCGSolver(); + /// Computes on-rank R and M_LH matrices std::pair ComputeSparseRAndM_LH(); + /// Abstract method to return the inverse of an on-rank lumped mass matrix virtual void LumpedMassInverse(Vector& ML_inv) const = 0; + /// Sets x_vdim_true, true dof values on a chosen vdim + /// given vector dof values x on the coarse (higher-order) mesh. virtual void GetHOTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const = 0; + /// Sets x_vdim_true, true dof values on a chosen vdim + /// given vector dof values x on the refined (low-order) mesh. virtual void GetLORTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const = 0; + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the coarse + /// (higher-order) mesh. virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const = 0; + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the refined + /// (low-order) mesh. virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const = 0; private: /// Computes sparsity pattern and initializes R matrix. Based on - /// BilinearForm::AllocMat() except maps between HO elements and LOR - /// elements. + /// BilinearForm::AllocMat() except maps between coarse HO elements and + /// refined LOR elements. SparseMatrix* AllocR(); }; + /** Implements 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) in serial. */ class SerialL2ProjectionH1Space : public L2ProjectionH1Space { private: @@ -337,12 +354,22 @@ protected: private: /// Computes inverse of a lumped mass matrix (stored as a vector) virtual void LumpedMassInverse(Vector& ML_inv) const; + /// Sets x_vdim_true, true dof values on a chosen vdim + /// given vector dof values x on the coarse (higher-order) mesh. virtual void GetHOTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const; + /// Sets x_vdim_true, true dof values on a chosen vdim + /// given vector dof values x on the refined (low-order) mesh. virtual void GetLORTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const; + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the coarse + /// (higher-order) mesh. virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const; + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the refined + /// (low-order) mesh. virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const; void GetDofsByVDim(const FiniteElementSpace& fes, @@ -353,6 +380,10 @@ protected: #ifdef MFEM_USE_MPI + + /** Implements 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) in parallel. */ class ParL2ProjectionH1Space : public L2ProjectionH1Space { private: @@ -371,12 +402,22 @@ protected: private: /// Computes inverse of a lumped mass matrix (stored as a vector) virtual void LumpedMassInverse(Vector& ML_inv) const; + /// Sets x_vdim_true, true dof values on a chosen vdim + /// given vector dof values x on the coarse (higher-order) mesh. virtual void GetHOTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const; + /// Sets x_vdim_true, true dof values on a chosen vdim + /// given vector dof values x on the refined (low-order) mesh. virtual void GetLORTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const; + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the coarse + /// (higher-order) mesh. virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const; + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the refined + /// (low-order) mesh. virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const; void GetTDofsByVDim(const ParFiniteElementSpace& pfes, From 35938b691ae8023d2fa8f78c0103da94850d9799 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 24 Mar 2023 12:58:02 -0700 Subject: [PATCH 07/34] astyle changes --- fem/transfer.cpp | 54 ++++++++++++++-------------- fem/transfer.hpp | 60 ++++++++++++++++---------------- miniapps/tools/lor-transferp.cpp | 21 ++++++----- 3 files changed, 71 insertions(+), 64 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index d0c1e52632..ec26b2f643 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -631,8 +631,8 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::InitializeCGSolver() SetAbsTol(1e-25); } -std::pair - L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() +std::pair +L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() { std::pair r_and_mlh; @@ -854,7 +854,7 @@ SparseMatrix* L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SerialL2ProjectionH1Space( const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_) -: L2ProjectionH1Space(fes_ho_, fes_lor_) + : L2ProjectionH1Space(fes_ho_, fes_lor_) { std::tie(R_sm, M_LH_sm) = ComputeSparseRAndM_LH(); RTxM_LH_sm = TransposeMult(*R_sm, *M_LH_sm); @@ -876,7 +876,8 @@ L2ProjectionGridTransfer::SerialL2ProjectionH1Space::~SerialL2ProjectionH1Space( delete RTxM_LH_sm; } -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::LumpedMassInverse(Vector& ML_inv) const +void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::LumpedMassInverse( + Vector& ML_inv) const { for (int i = 0; i < ML_inv.Size(); ++i) { @@ -901,9 +902,9 @@ void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::GetDofsByVDim( int vdim, const Vector& x, Vector& x_vdim) const { MFEM_ASSERT(x.Size() == fes.GetVSize(), - "Size of x Vector must match the number of vector DOFs."); + "Size of x Vector must match the number of vector DOFs."); MFEM_ASSERT(x_vdim.Size() == fes.GetNDofs(), - "Size of x_vdim Vector must match the number of DOFs."); + "Size of x_vdim Vector must match the number of DOFs."); Array vdim_dofs(fes.GetNDofs()); fes.GetVDofs(vdim, vdim_dofs); x.GetSubVector(vdim_dofs, x_vdim); @@ -926,9 +927,9 @@ void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SetFromDofsByVDim( int vdim, const Vector& y_vdim, Vector& y) const { MFEM_ASSERT(y_vdim.Size() == fes.GetNDofs(), - "Size of y_vdim Vector must match the number of DOFs."); + "Size of y_vdim Vector must match the number of DOFs."); MFEM_ASSERT(y.Size() == fes.GetVSize(), - "Size of y Vector must match the number of vector DOFs."); + "Size of y Vector must match the number of vector DOFs."); Array vdim_dofs(fes.GetNDofs()); fes.GetVDofs(vdim, vdim_dofs); y.SetSubVector(vdim_dofs, y_vdim); @@ -938,15 +939,15 @@ void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SetFromDofsByVDim( L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( const ParFiniteElementSpace& pfes_ho_, const ParFiniteElementSpace& pfes_lor_) -: L2ProjectionH1Space(pfes_ho_, pfes_lor_), - pfes_ho(pfes_ho_), - pfes_lor(pfes_lor_) + : L2ProjectionH1Space(pfes_ho_, pfes_lor_), + pfes_ho(pfes_ho_), + pfes_lor(pfes_lor_) { if (pfes_ho.GetVDim() != 1) { - pfes_ho_scalar = new ParFiniteElementSpace(pfes_ho.GetParMesh(), + pfes_ho_scalar = new ParFiniteElementSpace(pfes_ho.GetParMesh(), pfes_ho.FEColl(), 1); - pfes_lor_scalar = new ParFiniteElementSpace(pfes_lor.GetParMesh(), + pfes_lor_scalar = new ParFiniteElementSpace(pfes_lor.GetParMesh(), pfes_lor.FEColl(), 1); } else // vdim == 1 @@ -959,16 +960,16 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( std::tie(R_sm, M_LH_sm) = ComputeSparseRAndM_LH(); HypreParMatrix R_local = HypreParMatrix(pfes_ho_scalar->GetComm(), - pfes_lor_scalar->GlobalVSize(), pfes_ho_scalar->GlobalVSize(), - pfes_lor_scalar->GetDofOffsets(), pfes_ho_scalar->GetDofOffsets(), R_sm); + pfes_lor_scalar->GlobalVSize(), pfes_ho_scalar->GlobalVSize(), + pfes_lor_scalar->GetDofOffsets(), pfes_ho_scalar->GetDofOffsets(), R_sm); HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho_scalar->GetComm(), - pfes_lor_scalar->GlobalVSize(), pfes_ho_scalar->GlobalVSize(), - pfes_lor_scalar->GetDofOffsets(), pfes_ho_scalar->GetDofOffsets(), M_LH_sm); + pfes_lor_scalar->GlobalVSize(), pfes_ho_scalar->GlobalVSize(), + pfes_lor_scalar->GetDofOffsets(), pfes_ho_scalar->GetDofOffsets(), M_LH_sm); R_par = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(), &R_local, - pfes_ho_scalar->Dof_TrueDof_Matrix()); + pfes_ho_scalar->Dof_TrueDof_Matrix()); M_LH_par = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(), &M_LH_local, - pfes_ho_scalar->Dof_TrueDof_Matrix()); + pfes_ho_scalar->Dof_TrueDof_Matrix()); HypreParMatrix* R_T = R_par->Transpose(); RTxM_LH_par = ParMult(R_T, M_LH_par, true); @@ -1001,7 +1002,8 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::~ParL2ProjectionH1Space() delete M; } -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::LumpedMassInverse(Vector& ML_inv) const +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::LumpedMassInverse( + Vector& ML_inv) const { Vector ML_inv_true(pfes_lor_scalar->GetTrueVSize()); const Operator& P = *pfes_lor_scalar->GetProlongationMatrix(); @@ -1030,9 +1032,9 @@ void L2ProjectionGridTransfer::ParL2ProjectionH1Space::GetTDofsByVDim( int vdim, const Vector& x, Vector& x_vdim_true) const { MFEM_ASSERT(x.Size() == pfes.GetVSize(), - "Size of x Vector must match the number of vector DOFs."); + "Size of x Vector must match the number of vector DOFs."); MFEM_ASSERT(x_vdim_true.Size() == pfes_scalar.GetTrueVSize(), - "Size of x_vdim_true Vector must match the number of scalar true DOFs."); + "Size of x_vdim_true Vector must match the number of scalar true DOFs."); // transfer to vector of vdim dofs Array vdim_dofs(pfes.GetNDofs()); pfes.GetVDofs(vdim, vdim_dofs); @@ -1059,9 +1061,9 @@ void L2ProjectionGridTransfer::ParL2ProjectionH1Space::SetFromTDofsByVDim( int vdim, const Vector& y_vdim_true, Vector& y) const { MFEM_ASSERT(y_vdim_true.Size() == pfes_scalar.GetTrueVSize(), - "Size of y_vdim_true Vector must match the number of scalar true DOFs."); + "Size of y_vdim_true Vector must match the number of scalar true DOFs."); MFEM_ASSERT(y.Size() == pfes.GetVSize(), - "Size of y Vector must match the number of vector DOFs."); + "Size of y Vector must match the number of vector DOFs."); // prologate to vector of dofs on the vdim Vector y_vdim(pfes.GetNDofs()); pfes_scalar.GetProlongationMatrix()->Mult(y_vdim_true, y_vdim); @@ -1107,9 +1109,9 @@ void L2ProjectionGridTransfer::BuildF() else { #ifdef MFEM_USE_MPI - const mfem::ParFiniteElementSpace& dom_pfes = + const mfem::ParFiniteElementSpace& dom_pfes = static_cast(dom_fes); - const mfem::ParFiniteElementSpace& ran_pfes = + const mfem::ParFiniteElementSpace& ran_pfes = static_cast(ran_fes); F = new ParL2ProjectionH1Space(dom_pfes, ran_pfes); #endif diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 4f31c59c3b..42e7c7939d 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -314,21 +314,21 @@ protected: virtual void LumpedMassInverse(Vector& ML_inv) const = 0; /// Sets x_vdim_true, true dof values on a chosen vdim /// given vector dof values x on the coarse (higher-order) mesh. - virtual void GetHOTDofsByVDim(int vdim, const Vector& x, + virtual void GetHOTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const = 0; /// Sets x_vdim_true, true dof values on a chosen vdim /// given vector dof values x on the refined (low-order) mesh. - virtual void GetLORTDofsByVDim(int vdim, const Vector& x, + virtual void GetLORTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const = 0; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the coarse + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the coarse /// (higher-order) mesh. - virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const = 0; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the refined + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the refined /// (low-order) mesh. - virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const = 0; private: /// Computes sparsity pattern and initializes R matrix. Based on @@ -356,26 +356,26 @@ protected: virtual void LumpedMassInverse(Vector& ML_inv) const; /// Sets x_vdim_true, true dof values on a chosen vdim /// given vector dof values x on the coarse (higher-order) mesh. - virtual void GetHOTDofsByVDim(int vdim, const Vector& x, + virtual void GetHOTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const; /// Sets x_vdim_true, true dof values on a chosen vdim /// given vector dof values x on the refined (low-order) mesh. - virtual void GetLORTDofsByVDim(int vdim, const Vector& x, + virtual void GetLORTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the coarse + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the coarse /// (higher-order) mesh. - virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the refined + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the refined /// (low-order) mesh. - virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const; - void GetDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& x, Vector& x_vdim) const; - void SetFromDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& y_vdim, Vector& y) const; + void GetDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& x, Vector& x_vdim) const; + void SetFromDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& y_vdim, Vector& y) const; }; @@ -404,26 +404,26 @@ protected: virtual void LumpedMassInverse(Vector& ML_inv) const; /// Sets x_vdim_true, true dof values on a chosen vdim /// given vector dof values x on the coarse (higher-order) mesh. - virtual void GetHOTDofsByVDim(int vdim, const Vector& x, + virtual void GetHOTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const; /// Sets x_vdim_true, true dof values on a chosen vdim /// given vector dof values x on the refined (low-order) mesh. - virtual void GetLORTDofsByVDim(int vdim, const Vector& x, + virtual void GetLORTDofsByVDim(int vdim, const Vector& x, Vector& x_vdim_true) const; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the coarse + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the coarse /// (higher-order) mesh. - virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the refined + /// Sets the true dof values on a chosen vdim from + /// y_vdim_true onto y, vector dof values on the refined /// (low-order) mesh. - virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, + virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const; - void GetTDofsByVDim(const ParFiniteElementSpace& pfes, + void GetTDofsByVDim(const ParFiniteElementSpace& pfes, const ParFiniteElementSpace& pfes_scalar, int vdim, const Vector& x, Vector& x_vdim_true) const; - void SetFromTDofsByVDim(const ParFiniteElementSpace& pfes, + void SetFromTDofsByVDim(const ParFiniteElementSpace& pfes, const ParFiniteElementSpace& pfes_scalar, int vdim, const Vector& y_vdim_true, Vector& y) const; }; diff --git a/miniapps/tools/lor-transferp.cpp b/miniapps/tools/lor-transferp.cpp index abcd885cb3..ef21c109b9 100644 --- a/miniapps/tools/lor-transferp.cpp +++ b/miniapps/tools/lor-transferp.cpp @@ -229,13 +229,15 @@ int main(int argc, char *argv[]) M_rho.ParFESpace()->GetProlongationMatrix()->Mult(M_rho_true, M_rho); P.MultTranspose(M_rho, M_rho_lor); Vector M_rho_lor_true(M_rho_lor.ParFESpace()->GetTrueVSize()); - M_rho_lor.ParFESpace()->GetRestrictionOperator()->Mult(M_rho_lor, M_rho_lor_true); + M_rho_lor.ParFESpace()->GetRestrictionOperator()->Mult(M_rho_lor, + M_rho_lor_true); double local_ho_mass = M_rho_true.Sum(); double local_lor_mass = M_rho_lor_true.Sum(); double ho_mass; double lor_mass; MPI_Allreduce(&local_ho_mass, &ho_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&local_lor_mass, &lor_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&local_lor_mass, &lor_mass, 1, MPI_DOUBLE, MPI_SUM, + MPI_COMM_WORLD); if (Mpi::Root()) { cout << "HO -> LOR dual field: " << fabs(ho_mass - lor_mass) << endl << endl; @@ -285,7 +287,8 @@ int main(int argc, char *argv[]) rho_lor.GetTrueDofs(rho_lor_true); Vector M_rho_lor_true(M_rho_lor.ParFESpace()->GetTrueVSize()); M_lor_tdof->Mult(rho_lor_true, M_rho_lor_true); - M_rho_lor.ParFESpace()->GetProlongationMatrix()->Mult(M_rho_lor_true, M_rho_lor); + M_rho_lor.ParFESpace()->GetProlongationMatrix()->Mult(M_rho_lor_true, + M_rho_lor); R.MultTranspose(M_rho_lor, M_rho); Vector M_rho_true(M_rho.ParFESpace()->GetTrueVSize()); M_rho.ParFESpace()->GetRestrictionOperator()->Mult(M_rho, M_rho_true); @@ -294,7 +297,8 @@ int main(int argc, char *argv[]) double ho_mass; double lor_mass; MPI_Allreduce(&local_ho_mass, &ho_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&local_lor_mass, &lor_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&local_lor_mass, &lor_mass, 1, MPI_DOUBLE, MPI_SUM, + MPI_COMM_WORLD); if (Mpi::Root()) { cout << "LOR -> HO dual field: " << fabs(ho_mass - lor_mass) << '\n'; @@ -339,7 +343,8 @@ void visualize(VisItDataCollection &dc, string prefix, int x, int y) int visport = 19916; socketstream sol_sockL2(vishost, visport); - sol_sockL2 << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << "\n"; + 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 @@ -360,7 +365,7 @@ double compute_mass(ParFiniteElementSpace *L2, double massL2, Vector rhoone(L2->GetTrueVSize()); rhoone = 1.0; - + Vector Mdiag(L2->GetTrueVSize()); pML2->Mult(rhoone, Mdiag); delete pML2; @@ -373,8 +378,8 @@ double compute_mass(ParFiniteElementSpace *L2, double massL2, cout << space << " " << prefix << " mass = " << newmass; if (massL2 >= 0) { - cout.precision(4); - cout << " (" << fabs(newmass-massL2)*100/massL2 << "%)"; + cout.precision(4); + cout << " (" << fabs(newmass-massL2)*100/massL2 << "%)"; } cout << endl; } From b81b023ea55741c0dec1158f7b890eed112cd1e2 Mon Sep 17 00:00:00 2001 From: EB Chin Date: Fri, 24 Mar 2023 13:35:57 -0700 Subject: [PATCH 08/34] simplify parallel get/set dofs --- fem/transfer.cpp | 64 +++++++++++++++++++++--------------------------- fem/transfer.hpp | 8 +++--- 2 files changed, 32 insertions(+), 40 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index d0c1e52632..113766dcc8 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -763,6 +763,32 @@ std::pair return r_and_mlh; } +void L2ProjectionGridTransfer::L2ProjectionH1Space::GetDofsByVDim( + const FiniteElementSpace& fes, + int vdim, const Vector& x, Vector& x_vdim) const +{ + MFEM_ASSERT(x.Size() == fes.GetVSize(), + "Size of x Vector must match the number of vector DOFs."); + MFEM_ASSERT(x_vdim.Size() == fes.GetNDofs(), + "Size of x_vdim Vector must match the number of DOFs."); + Array vdim_dofs(fes.GetNDofs()); + fes.GetVDofs(vdim, vdim_dofs); + x.GetSubVector(vdim_dofs, x_vdim); +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromDofsByVDim( + const FiniteElementSpace& fes, + int vdim, const Vector& y_vdim, Vector& y) const +{ + MFEM_ASSERT(y_vdim.Size() == fes.GetNDofs(), + "Size of y_vdim Vector must match the number of DOFs."); + MFEM_ASSERT(y.Size() == fes.GetVSize(), + "Size of y Vector must match the number of vector DOFs."); + Array vdim_dofs(fes.GetNDofs()); + fes.GetVDofs(vdim, vdim_dofs); + y.SetSubVector(vdim_dofs, y_vdim); +} + SparseMatrix* L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() { const Table& elem_dof_ho = fes_ho.GetElementToDofTable(); @@ -896,19 +922,6 @@ void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::GetLORTDofsByVDim( GetDofsByVDim(fes_lor, vdim, x, x_vdim_true); } -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::GetDofsByVDim( - const FiniteElementSpace& fes, - int vdim, const Vector& x, Vector& x_vdim) const -{ - MFEM_ASSERT(x.Size() == fes.GetVSize(), - "Size of x Vector must match the number of vector DOFs."); - MFEM_ASSERT(x_vdim.Size() == fes.GetNDofs(), - "Size of x_vdim Vector must match the number of DOFs."); - Array vdim_dofs(fes.GetNDofs()); - fes.GetVDofs(vdim, vdim_dofs); - x.GetSubVector(vdim_dofs, x_vdim); -} - void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SetHOFromTDofsByVDim( int vdim, const Vector& y_vdim_true, Vector& y) const { @@ -921,19 +934,6 @@ void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SetLORFromTDofsByVDim( SetFromDofsByVDim(fes_lor, vdim, y_vdim_true, y); } -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SetFromDofsByVDim( - const FiniteElementSpace& fes, - int vdim, const Vector& y_vdim, Vector& y) const -{ - MFEM_ASSERT(y_vdim.Size() == fes.GetNDofs(), - "Size of y_vdim Vector must match the number of DOFs."); - MFEM_ASSERT(y.Size() == fes.GetVSize(), - "Size of y Vector must match the number of vector DOFs."); - Array vdim_dofs(fes.GetNDofs()); - fes.GetVDofs(vdim, vdim_dofs); - y.SetSubVector(vdim_dofs, y_vdim); -} - #ifdef MFEM_USE_MPI L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( @@ -1029,15 +1029,11 @@ void L2ProjectionGridTransfer::ParL2ProjectionH1Space::GetTDofsByVDim( const ParFiniteElementSpace& pfes, const ParFiniteElementSpace& pfes_scalar, int vdim, const Vector& x, Vector& x_vdim_true) const { - MFEM_ASSERT(x.Size() == pfes.GetVSize(), - "Size of x Vector must match the number of vector DOFs."); MFEM_ASSERT(x_vdim_true.Size() == pfes_scalar.GetTrueVSize(), "Size of x_vdim_true Vector must match the number of scalar true DOFs."); // transfer to vector of vdim dofs - Array vdim_dofs(pfes.GetNDofs()); - pfes.GetVDofs(vdim, vdim_dofs); Vector x_vdim(pfes.GetNDofs()); - x.GetSubVector(vdim_dofs, x_vdim); + GetDofsByVDim(pfes, vdim, x, x_vdim); // transfer to vector of true dofs on the vdim pfes_scalar.GetRestrictionOperator()->Mult(x_vdim, x_vdim_true); } @@ -1060,15 +1056,11 @@ void L2ProjectionGridTransfer::ParL2ProjectionH1Space::SetFromTDofsByVDim( { MFEM_ASSERT(y_vdim_true.Size() == pfes_scalar.GetTrueVSize(), "Size of y_vdim_true Vector must match the number of scalar true DOFs."); - MFEM_ASSERT(y.Size() == pfes.GetVSize(), - "Size of y Vector must match the number of vector DOFs."); // prologate to vector of dofs on the vdim Vector y_vdim(pfes.GetNDofs()); pfes_scalar.GetProlongationMatrix()->Mult(y_vdim_true, y_vdim); // transfer to global vector - Array vdim_dofs(pfes.GetNDofs()); - pfes.GetVDofs(vdim, vdim_dofs); - y.SetSubVector(vdim_dofs, y_vdim); + SetFromDofsByVDim(pfes, vdim, y_vdim, y); } #endif diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 4f31c59c3b..1d0e4a3dac 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -310,6 +310,10 @@ protected: void InitializeCGSolver(); /// Computes on-rank R and M_LH matrices std::pair ComputeSparseRAndM_LH(); + void GetDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& x, Vector& x_vdim) const; + void SetFromDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& y_vdim, Vector& y) const; /// Abstract method to return the inverse of an on-rank lumped mass matrix virtual void LumpedMassInverse(Vector& ML_inv) const = 0; /// Sets x_vdim_true, true dof values on a chosen vdim @@ -372,10 +376,6 @@ protected: /// (low-order) mesh. virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, Vector& y) const; - void GetDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& x, Vector& x_vdim) const; - void SetFromDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& y_vdim, Vector& y) const; }; From 102d4ad3f8ad4ec9b2662d09ba2c2da79619c4aa Mon Sep 17 00:00:00 2001 From: EB Chin Date: Fri, 24 Mar 2023 13:44:44 -0700 Subject: [PATCH 09/34] fix output --- miniapps/tools/lor-transferp.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/miniapps/tools/lor-transferp.cpp b/miniapps/tools/lor-transferp.cpp index ef21c109b9..0de4baf5fb 100644 --- a/miniapps/tools/lor-transferp.cpp +++ b/miniapps/tools/lor-transferp.cpp @@ -131,7 +131,10 @@ int main(int argc, char *argv[]) if (lorder == 0) { lorder = 1; - cerr << "Switching the H1 LOR space order from 0 to 1\n"; + 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); From 36e3062378f44b25b70de427a2d8aea2da912b17 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 24 Mar 2023 13:47:06 -0700 Subject: [PATCH 10/34] make style --- fem/transfer.cpp | 8 ++++---- fem/transfer.hpp | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 5905afed79..2d937db319 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -767,9 +767,9 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::GetDofsByVDim( int vdim, const Vector& x, Vector& x_vdim) const { MFEM_ASSERT(x.Size() == fes.GetVSize(), - "Size of x Vector must match the number of vector DOFs."); + "Size of x Vector must match the number of vector DOFs."); MFEM_ASSERT(x_vdim.Size() == fes.GetNDofs(), - "Size of x_vdim Vector must match the number of DOFs."); + "Size of x_vdim Vector must match the number of DOFs."); Array vdim_dofs(fes.GetNDofs()); fes.GetVDofs(vdim, vdim_dofs); x.GetSubVector(vdim_dofs, x_vdim); @@ -780,9 +780,9 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromDofsByVDim( int vdim, const Vector& y_vdim, Vector& y) const { MFEM_ASSERT(y_vdim.Size() == fes.GetNDofs(), - "Size of y_vdim Vector must match the number of DOFs."); + "Size of y_vdim Vector must match the number of DOFs."); MFEM_ASSERT(y.Size() == fes.GetVSize(), - "Size of y Vector must match the number of vector DOFs."); + "Size of y Vector must match the number of vector DOFs."); Array vdim_dofs(fes.GetNDofs()); fes.GetVDofs(vdim, vdim_dofs); y.SetSubVector(vdim_dofs, y_vdim); diff --git a/fem/transfer.hpp b/fem/transfer.hpp index a3fdf382d0..c6342aee0d 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -310,10 +310,10 @@ protected: void InitializeCGSolver(); /// Computes on-rank R and M_LH matrices std::pair ComputeSparseRAndM_LH(); - void GetDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& x, Vector& x_vdim) const; - void SetFromDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& y_vdim, Vector& y) const; + void GetDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& x, Vector& x_vdim) const; + void SetFromDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& y_vdim, Vector& y) const; /// Abstract method to return the inverse of an on-rank lumped mass matrix virtual void LumpedMassInverse(Vector& ML_inv) const = 0; /// Sets x_vdim_true, true dof values on a chosen vdim From e21949b16baf54f17198d06d4755eb4c1ab3af49 Mon Sep 17 00:00:00 2001 From: EB Chin Date: Fri, 24 Mar 2023 20:58:00 -0700 Subject: [PATCH 11/34] update copyright date --- miniapps/tools/lor-transferp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/miniapps/tools/lor-transferp.cpp b/miniapps/tools/lor-transferp.cpp index 0de4baf5fb..cabe13ddf1 100644 --- a/miniapps/tools/lor-transferp.cpp +++ b/miniapps/tools/lor-transferp.cpp @@ -1,4 +1,4 @@ -// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced +// 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. // From 2aa15a7ef0b0463197d4d024d36a9da678f441d6 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Thu, 13 Jul 2023 17:54:33 -0700 Subject: [PATCH 12/34] address review comments --- .gitignore | 1 + fem/transfer.cpp | 414 +++++++++--------- fem/transfer.hpp | 147 ++----- miniapps/tools/CMakeLists.txt | 4 +- miniapps/tools/lor-transfer.cpp | 8 +- miniapps/tools/makefile | 6 +- .../{lor-transferp.cpp => plor-transfer.cpp} | 72 ++- 7 files changed, 266 insertions(+), 386 deletions(-) rename miniapps/tools/{lor-transferp.cpp => plor-transfer.cpp} (87%) diff --git a/.gitignore b/.gitignore index 78f76d12e6..db9034b959 100644 --- a/.gitignore +++ b/.gitignore @@ -286,6 +286,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/fem/transfer.cpp b/fem/transfer.cpp index 5c355c4926..677bfbc120 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -544,96 +544,136 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::ProlongateTranspose( L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_) - : L2Projection(fes_ho_, fes_lor_) -{} - -L2ProjectionGridTransfer::L2ProjectionH1Space::~L2ProjectionH1Space() -{} + : L2Projection(fes_ho_, fes_lor_), +#ifdef MFEM_USE_MPI + pcg { dynamic_cast(&fes_ho) ? + CGSolver(MPI_COMM_WORLD) : CGSolver() } +#else + pcg() +#endif +{ + std::tie(R, M_LH) = ComputeSparseRAndM_LH(); + RTxM_LH = std::unique_ptr(TransposeMult(static_cast(*R), + static_cast(*M_LH))); + + // Basic PCG solver setup + pcg.SetPrintLevel(0); + pcg.SetMaxIter(1000); + // initial values for relative and absolute tolerance + SetRelTol(1e-13); + SetAbsTol(1e-13); + precon = std::unique_ptr(new DSmoother(static_cast(*RTxM_LH))); + 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()); + + GetTDofs(fes_ho, x, X); + for (int d = 0; d < fes_ho.GetVDim(); ++d) { - GetHOTDofsByVDim(d, x, X_dim); + GetTDofsByVDim(fes_ho, d, X, X_dim); R->Mult(X_dim, Y_dim); - SetLORFromTDofsByVDim(d, Y_dim, y); + SetTDofsByVDim(fes_lor, d, Y_dim, Y); } + + SetTDofs(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()); + + GetTDofs(fes_lor, x, X); + for (int d = 0; d < fes_ho.GetVDim(); ++d) { - GetLORTDofsByVDim(d, x, X_dim); + GetTDofsByVDim(fes_lor, d, X, X_dim); R->MultTranspose(X_dim, Y_dim); - SetHOFromTDofsByVDim(d, Y_dim, y); + SetTDofsByVDim(fes_ho, d, Y_dim, Y); } + + SetTDofs(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 Y_dim(pcg.Height()); Vector Xbar(pcg.Width()); + + Vector Y_dim(pcg.Height()); + Vector Y(fes_ho.GetTrueVSize()); + + GetTDofs(fes_lor, x, X); + for (int d = 0; d < fes_ho.GetVDim(); ++d) { - GetLORTDofsByVDim(d, x, X_dim); - // Compute y = P x = (R^T M_LH)^(-1) M_LH^T x = (R^T M_LH)^(-1) xbar + GetTDofsByVDim(fes_lor, d, X, 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); - SetHOFromTDofsByVDim(d, Y_dim, y); + SetTDofsByVDim(fes_ho, d, Y_dim, Y); } + + SetTDofs(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 Y_dim(M_LH->Height()); Vector Xbar(pcg.Height()); + + Vector Y_dim(M_LH->Height()); + Vector Y(fes_lor.GetTrueVSize()); + + GetTDofs(fes_ho, x, X); + for (int d = 0; d < fes_ho.GetVDim(); ++d) { - GetHOTDofsByVDim(d, x, X_dim); - // Compute y = P^T x = M_LH (R^T M_LH)^(-1) x = M_LH xbar + GetTDofsByVDim(fes_ho, d, X, 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); - SetLORFromTDofsByVDim(d, Y_dim, y); + SetTDofsByVDim(fes_lor, d, Y_dim, y); } + + SetTDofs(fes_lor, Y, y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_) { - pcg.SetRelTol(sqrt(p_rtol_)); + pcg.SetRelTol(p_rtol_); } void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_) { - pcg.SetAbsTol(sqrt(p_atol_)); + pcg.SetAbsTol(p_atol_); } -void L2ProjectionGridTransfer::L2ProjectionH1Space::InitializeCGSolver() -{ - // Basic PCG solver setup - pcg.SetPrintLevel(0); - pcg.SetMaxIter(1000); - // initial values for relative and absolute tolerance - SetRelTol(1e-25); - SetAbsTol(1e-25); -} - -std::pair +std::pair, std::unique_ptr> L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() { - std::pair r_and_mlh; + std::pair, std::unique_ptr> + r_and_mlh; Mesh* mesh_ho = fes_ho.GetMesh(); Mesh* mesh_lor = fes_lor.GetMesh(); @@ -707,9 +747,9 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() // 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) - r_and_mlh.second = new SparseMatrix( + r_and_mlh.second = std::unique_ptr(new SparseMatrix( r_and_mlh.first->GetI(), r_and_mlh.first->GetJ(), NULL, - r_and_mlh.first->Height(), r_and_mlh.first->Width(), false, true, true); + r_and_mlh.first->Height(), r_and_mlh.first->Width(), false, true, true)); IntegrationPointTransformation ip_tr; IsoparametricTransformation& emb_tr = ip_tr.Transf; @@ -762,33 +802,62 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() return r_and_mlh; } -void L2ProjectionGridTransfer::L2ProjectionH1Space::GetDofsByVDim( - const FiniteElementSpace& fes, - int vdim, const Vector& x, Vector& x_vdim) const +void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofs( + const FiniteElementSpace& fes, const Vector& x, Vector& X) const { - MFEM_ASSERT(x.Size() == fes.GetVSize(), - "Size of x Vector must match the number of vector DOFs."); - MFEM_ASSERT(x_vdim.Size() == fes.GetNDofs(), - "Size of x_vdim Vector must match the number of DOFs."); - Array vdim_dofs(fes.GetNDofs()); - fes.GetVDofs(vdim, vdim_dofs); - x.GetSubVector(vdim_dofs, x_vdim); + const Operator* res = fes.GetRestrictionOperator(); + if (res) + { + res->Mult(x, X); + } + else + { + X.NewDataAndSize(x.GetData(), x.Size()); + } } -void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromDofsByVDim( - const FiniteElementSpace& fes, - int vdim, const Vector& y_vdim, Vector& y) const +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofs( + const FiniteElementSpace& fes, const Vector& X, Vector& x) const { - MFEM_ASSERT(y_vdim.Size() == fes.GetNDofs(), - "Size of y_vdim Vector must match the number of DOFs."); - MFEM_ASSERT(y.Size() == fes.GetVSize(), - "Size of y Vector must match the number of vector DOFs."); - Array vdim_dofs(fes.GetNDofs()); - fes.GetVDofs(vdim, vdim_dofs); - y.SetSubVector(vdim_dofs, y_vdim); + const Operator* P = fes.GetProlongationMatrix(); + if (P) + { + P->Mult(X, x); + } + else + { + x.NewDataAndSize(X.GetData(), X.Size()); + } } -SparseMatrix* L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() +void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofsByVDim( + const FiniteElementSpace& fes, int vdim, const Vector& X, Vector& X_vdim) const +{ + Array x_vdofs_list(fes.GetNDofs()); + + fes.GetVDofs(vdim, x_vdofs_list); + X.GetSubVector(x_vdofs_list, X_vdim); +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofsByVDim( + const FiniteElementSpace& fes, int vdim, const Vector& X_vdim, Vector& X) const +{ + Array x_vdofs_list(fes.GetNDofs()); + + fes.GetVDofs(vdim, x_vdofs_list); + X.SetSubVector(x_vdofs_list, X_vdim); +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::LumpedMassInverse( + Vector& ML_inv) const +{ + for (int i = 0; i < ML_inv.Size(); ++i) + { + ML_inv[i] = 1.0 / ML_inv[i]; + } +} + +std::unique_ptr L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() { const Table& elem_dof_ho = fes_ho.GetElementToDofTable(); const Table& elem_dof_lor = fes_lor.GetElementToDofTable(); @@ -868,8 +937,8 @@ SparseMatrix* L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() dof_lor_dof_ho.SortRows(); double* data = Memory(dof_dofI[ndof_lor]); - SparseMatrix* R_local = new SparseMatrix( - dof_dofI, dof_dofJ, data, ndof_lor, ndof_ho, true, true, true); + 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(); @@ -877,63 +946,6 @@ SparseMatrix* L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() return R_local; } -L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SerialL2ProjectionH1Space( - const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_) - : L2ProjectionH1Space(fes_ho_, fes_lor_) -{ - std::tie(R_sm, M_LH_sm) = ComputeSparseRAndM_LH(); - RTxM_LH_sm = TransposeMult(*R_sm, *M_LH_sm); - - R = R_sm; - M_LH = M_LH_sm; - RTxM_LH = RTxM_LH_sm; - - InitializeCGSolver(); - Ds = DSmoother(*RTxM_LH_sm); - pcg.SetPreconditioner(Ds); - pcg.SetOperator(*RTxM_LH); -} - -L2ProjectionGridTransfer::SerialL2ProjectionH1Space::~SerialL2ProjectionH1Space() -{ - delete R_sm; - delete M_LH_sm; - delete RTxM_LH_sm; -} - -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::LumpedMassInverse( - Vector& ML_inv) const -{ - for (int i = 0; i < ML_inv.Size(); ++i) - { - ML_inv[i] = 1.0 / ML_inv[i]; - } -} - -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::GetHOTDofsByVDim( - int vdim, const Vector& x, Vector& x_vdim_true) const -{ - GetDofsByVDim(fes_ho, vdim, x, x_vdim_true); -} - -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::GetLORTDofsByVDim( - int vdim, const Vector& x, Vector& x_vdim_true) const -{ - GetDofsByVDim(fes_lor, vdim, x, x_vdim_true); -} - -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SetHOFromTDofsByVDim( - int vdim, const Vector& y_vdim_true, Vector& y) const -{ - SetFromDofsByVDim(fes_ho, vdim, y_vdim_true, y); -} - -void L2ProjectionGridTransfer::SerialL2ProjectionH1Space::SetLORFromTDofsByVDim( - int vdim, const Vector& y_vdim_true, Vector& y) const -{ - SetFromDofsByVDim(fes_lor, vdim, y_vdim_true, y); -} - #ifdef MFEM_USE_MPI L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( @@ -942,126 +954,94 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( pfes_ho(pfes_ho_), pfes_lor(pfes_lor_) { - if (pfes_ho.GetVDim() != 1) - { - pfes_ho_scalar = new ParFiniteElementSpace(pfes_ho.GetParMesh(), - pfes_ho.FEColl(), 1); - pfes_lor_scalar = new ParFiniteElementSpace(pfes_lor.GetParMesh(), - pfes_lor.FEColl(), 1); - } - else // vdim == 1 - { - pfes_ho_scalar = &pfes_ho; - pfes_lor_scalar = &pfes_lor; - } - SparseMatrix* R_sm; - SparseMatrix* M_LH_sm; - std::tie(R_sm, M_LH_sm) = 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_scalar->GetComm(), - pfes_lor_scalar->GlobalVSize(), pfes_ho_scalar->GlobalVSize(), - pfes_lor_scalar->GetDofOffsets(), pfes_ho_scalar->GetDofOffsets(), R_sm); - HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho_scalar->GetComm(), - pfes_lor_scalar->GlobalVSize(), pfes_ho_scalar->GlobalVSize(), - pfes_lor_scalar->GetDofOffsets(), pfes_ho_scalar->GetDofOffsets(), M_LH_sm); + 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())); - R_par = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(), &R_local, - pfes_ho_scalar->Dof_TrueDof_Matrix()); - M_LH_par = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(), &M_LH_local, - pfes_ho_scalar->Dof_TrueDof_Matrix()); - HypreParMatrix* R_T = R_par->Transpose(); - RTxM_LH_par = ParMult(R_T, M_LH_par, true); + R = std::unique_ptr(RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), + &R_local, pfes_ho_scalar.Dof_TrueDof_Matrix())); + M_LH = std::unique_ptr(RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), + &M_LH_local, pfes_ho_scalar.Dof_TrueDof_Matrix())); + std::unique_ptr R_T = + std::unique_ptr(static_cast(*R).Transpose()); + RTxM_LH = std::unique_ptr(ParMult( + R_T.get(), static_cast(M_LH.get()), true)); - R = R_par; - M_LH = M_LH_par; - RTxM_LH = RTxM_LH_par; - - pcg = CGSolver(MPI_COMM_WORLD); - InitializeCGSolver(); - M = new HypreBoomerAMG(*RTxM_LH_par); - M->SetPrintLevel(0); - pcg.SetPreconditioner(*M); - pcg.SetOperator(*RTxM_LH_par); - - delete R_sm; - delete M_LH_sm; - delete R_T; + precon = std::unique_ptr( + new HypreBoomerAMG(static_cast(*RTxM_LH))); + static_cast(*precon).SetPrintLevel(0); + pcg.SetPreconditioner(*precon); + pcg.SetOperator(*RTxM_LH); } -L2ProjectionGridTransfer::ParL2ProjectionH1Space::~ParL2ProjectionH1Space() +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::GetTDofsByVDim( + const FiniteElementSpace& fes, int vdim, const Vector& X, Vector& X_vdim) const { - if (pfes_ho.GetVDim() != 1) - { - delete pfes_ho_scalar; - delete pfes_lor_scalar; - } - delete R_par; - delete M_LH_par; - delete RTxM_LH_par; - delete M; + Array x_vdofs_list(fes.GetNDofs()); + Array x_vdofs_marker(fes.GetVSize()); + Array X_vdofs_marker(fes.GetTrueVSize()); + Array X_vdofs_list; + + fes.GetVDofs(vdim, x_vdofs_list); + FiniteElementSpace::ListToMarker(x_vdofs_list, fes.GetVSize(), x_vdofs_marker); + fes.GetRestrictionMatrix()->BooleanMult(x_vdofs_marker, X_vdofs_marker); + FiniteElementSpace::MarkerToList(X_vdofs_marker, X_vdofs_list); + X.GetSubVector(X_vdofs_list, X_vdim); +} + +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::SetTDofsByVDim( + const FiniteElementSpace& fes, int vdim, const Vector& X_vdim, Vector& X) const +{ + const ParFiniteElementSpace* pfes = dynamic_cast(&fes); + MFEM_ASSERT(pfes != nullptr, "fes must be a parallel space"); + + Array x_vdofs_list(fes.GetNDofs()); + Array x_vdofs_marker(fes.GetVSize()); + Array X_vdofs_marker(fes.GetTrueVSize()); + Array X_vdofs_list; + + fes.GetVDofs(vdim, x_vdofs_list); + FiniteElementSpace::ListToMarker(x_vdofs_list, fes.GetVSize(), x_vdofs_marker); + pfes->Dof_TrueDof_Matrix()->BooleanMult(1.0, x_vdofs_marker, 0.0, X_vdofs_marker); + FiniteElementSpace::MarkerToList(X_vdofs_marker, X_vdofs_list); + X.SetSubVector(X_vdofs_list, X_vdim); } void L2ProjectionGridTransfer::ParL2ProjectionH1Space::LumpedMassInverse( Vector& ML_inv) const { - Vector ML_inv_true(pfes_lor_scalar->GetTrueVSize()); - const Operator& P = *pfes_lor_scalar->GetProlongationMatrix(); - P.MultTranspose(ML_inv, ML_inv_true); + Vector ML_inv_full(pfes_lor.GetVSize()); + // set ML_inv on dofs for vdim = 0 + Array vdofs_list(pfes_lor.GetNDofs()); + pfes_lor.GetVDofs(0, vdofs_list); + ML_inv_full.SetSubVector(vdofs_list, ML_inv); + + Vector ML_inv_true(pfes_lor.GetTrueVSize()); + const Operator& P = *pfes_lor.GetProlongationMatrix(); + P.MultTranspose(ML_inv_full, ML_inv_true); + for (int i = 0; i < ML_inv_true.Size(); ++i) { ML_inv_true[i] = 1.0 / ML_inv_true[i]; } - P.Mult(ML_inv_true, ML_inv); -} -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::GetHOTDofsByVDim( - int vdim, const Vector& x, Vector& x_vdim_true) const -{ - GetTDofsByVDim(pfes_ho, *pfes_ho_scalar, vdim, x, x_vdim_true); -} + P.Mult(ML_inv_true, ML_inv_full); -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::GetLORTDofsByVDim( - int vdim, const Vector& x, Vector& x_vdim_true) const -{ - GetTDofsByVDim(pfes_lor, *pfes_lor_scalar, vdim, x, x_vdim_true); -} - -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::GetTDofsByVDim( - const ParFiniteElementSpace& pfes, const ParFiniteElementSpace& pfes_scalar, - int vdim, const Vector& x, Vector& x_vdim_true) const -{ - MFEM_ASSERT(x_vdim_true.Size() == pfes_scalar.GetTrueVSize(), - "Size of x_vdim_true Vector must match the number of scalar true DOFs."); - // transfer to vector of vdim dofs - Vector x_vdim(pfes.GetNDofs()); - GetDofsByVDim(pfes, vdim, x, x_vdim); - // transfer to vector of true dofs on the vdim - pfes_scalar.GetRestrictionOperator()->Mult(x_vdim, x_vdim_true); -} - -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::SetHOFromTDofsByVDim( - int vdim, const Vector& y_vdim_true, Vector& y) const -{ - SetFromTDofsByVDim(pfes_ho, *pfes_ho_scalar, vdim, y_vdim_true, y); -} - -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::SetLORFromTDofsByVDim( - int vdim, const Vector& y_vdim_true, Vector& y) const -{ - SetFromTDofsByVDim(pfes_lor, *pfes_lor_scalar, vdim, y_vdim_true, y); -} - -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::SetFromTDofsByVDim( - const ParFiniteElementSpace& pfes, const ParFiniteElementSpace& pfes_scalar, - int vdim, const Vector& y_vdim_true, Vector& y) const -{ - MFEM_ASSERT(y_vdim_true.Size() == pfes_scalar.GetTrueVSize(), - "Size of y_vdim_true Vector must match the number of scalar true DOFs."); - // prologate to vector of dofs on the vdim - Vector y_vdim(pfes.GetNDofs()); - pfes_scalar.GetProlongationMatrix()->Mult(y_vdim_true, y_vdim); - // transfer to global vector - SetFromDofsByVDim(pfes, vdim, y_vdim, y); + ML_inv_full.GetSubVector(vdofs_list, ML_inv); } #endif @@ -1095,7 +1075,7 @@ void L2ProjectionGridTransfer::BuildF() { if (!Parallel()) { - F = new SerialL2ProjectionH1Space(dom_fes, ran_fes); + F = new L2ProjectionH1Space(dom_fes, ran_fes); } else { diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 7a60359d4e..f6d7ef3bb2 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -258,20 +258,10 @@ protected: on a refined mesh (LOR). */ class L2ProjectionH1Space : public L2Projection { - protected: - // 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. - Operator* R; - // Used to compute P = (RTxM_LH)^(-1) M_LH^T - Operator* M_LH; - Operator* RTxM_LH; - CGSolver pcg; - + public: L2ProjectionH1Space(const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_); - public: - virtual ~L2ProjectionH1Space(); + virtual ~L2ProjectionH1Space() = default; /// 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 @@ -306,76 +296,34 @@ protected: virtual void SetRelTol(double p_rtol_); virtual void SetAbsTol(double p_atol_); protected: - /// Sets default shared options for the CGSolver - void InitializeCGSolver(); /// Computes on-rank R and M_LH matrices - std::pair ComputeSparseRAndM_LH(); - void GetDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& x, Vector& x_vdim) const; - void SetFromDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& y_vdim, Vector& y) const; - /// Abstract method to return the inverse of an on-rank lumped mass matrix - virtual void LumpedMassInverse(Vector& ML_inv) const = 0; - /// Sets x_vdim_true, true dof values on a chosen vdim - /// given vector dof values x on the coarse (higher-order) mesh. - virtual void GetHOTDofsByVDim(int vdim, const Vector& x, - Vector& x_vdim_true) const = 0; - /// Sets x_vdim_true, true dof values on a chosen vdim - /// given vector dof values x on the refined (low-order) mesh. - virtual void GetLORTDofsByVDim(int vdim, const Vector& x, - Vector& x_vdim_true) const = 0; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the coarse - /// (higher-order) mesh. - virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, - Vector& y) const = 0; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the refined - /// (low-order) mesh. - virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, - Vector& y) const = 0; + std::pair, std::unique_ptr> + ComputeSparseRAndM_LH(); + virtual void GetTDofs(const FiniteElementSpace& fes, + const Vector& x, Vector& X) const; + virtual void SetTDofs(const FiniteElementSpace& fes, + const Vector& X, Vector& x) const; + virtual void GetTDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& X, Vector& X_vdim) const; + virtual void SetTDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& X_vdim, Vector& X) const; + /// Returns the inverse of an on-rank lumped mass matrix + virtual void LumpedMassInverse(Vector& ML_inv) const; + + 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; private: /// Computes sparsity pattern and initializes R matrix. Based on /// BilinearForm::AllocMat() except maps between coarse HO elements and /// refined LOR elements. - SparseMatrix* AllocR(); - }; - - /** Implements 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) in serial. */ - class SerialL2ProjectionH1Space : public L2ProjectionH1Space - { - private: - SparseMatrix* R_sm; - SparseMatrix* M_LH_sm; - SparseMatrix* RTxM_LH_sm; - DSmoother Ds; - public: - SerialL2ProjectionH1Space(const FiniteElementSpace& fes_ho_, - const FiniteElementSpace& fes_lor_); - virtual ~SerialL2ProjectionH1Space(); - private: - /// Computes inverse of a lumped mass matrix (stored as a vector) - virtual void LumpedMassInverse(Vector& ML_inv) const; - /// Sets x_vdim_true, true dof values on a chosen vdim - /// given vector dof values x on the coarse (higher-order) mesh. - virtual void GetHOTDofsByVDim(int vdim, const Vector& x, - Vector& x_vdim_true) const; - /// Sets x_vdim_true, true dof values on a chosen vdim - /// given vector dof values x on the refined (low-order) mesh. - virtual void GetLORTDofsByVDim(int vdim, const Vector& x, - Vector& x_vdim_true) const; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the coarse - /// (higher-order) mesh. - virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, - Vector& y) const; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the refined - /// (low-order) mesh. - virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, - Vector& y) const; + std::unique_ptr AllocR(); }; @@ -386,47 +334,22 @@ protected: on a refined mesh (LOR) in parallel. */ class ParL2ProjectionH1Space : public L2ProjectionH1Space { - private: - const ParFiniteElementSpace& pfes_ho; - const ParFiniteElementSpace& pfes_lor; - const ParFiniteElementSpace* pfes_ho_scalar; - const ParFiniteElementSpace* pfes_lor_scalar; - HypreParMatrix* R_par; - HypreParMatrix* M_LH_par; - HypreParMatrix* RTxM_LH_par; - HypreBoomerAMG* M; public: ParL2ProjectionH1Space(const ParFiniteElementSpace& pfes_ho_, const ParFiniteElementSpace& pfes_lor_); - virtual ~ParL2ProjectionH1Space(); + virtual ~ParL2ProjectionH1Space() = default; private: + void GetTDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& x, Vector& x_vdim) const override; + void SetTDofsByVDim(const FiniteElementSpace& fes, + int vdim, const Vector& x_vdim, Vector& x) const override; /// Computes inverse of a lumped mass matrix (stored as a vector) - virtual void LumpedMassInverse(Vector& ML_inv) const; - /// Sets x_vdim_true, true dof values on a chosen vdim - /// given vector dof values x on the coarse (higher-order) mesh. - virtual void GetHOTDofsByVDim(int vdim, const Vector& x, - Vector& x_vdim_true) const; - /// Sets x_vdim_true, true dof values on a chosen vdim - /// given vector dof values x on the refined (low-order) mesh. - virtual void GetLORTDofsByVDim(int vdim, const Vector& x, - Vector& x_vdim_true) const; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the coarse - /// (higher-order) mesh. - virtual void SetHOFromTDofsByVDim(int vdim, const Vector& y_vdim_true, - Vector& y) const; - /// Sets the true dof values on a chosen vdim from - /// y_vdim_true onto y, vector dof values on the refined - /// (low-order) mesh. - virtual void SetLORFromTDofsByVDim(int vdim, const Vector& y_vdim_true, - Vector& y) const; - void GetTDofsByVDim(const ParFiniteElementSpace& pfes, - const ParFiniteElementSpace& pfes_scalar, - int vdim, const Vector& x, Vector& x_vdim_true) const; - void SetFromTDofsByVDim(const ParFiniteElementSpace& pfes, - const ParFiniteElementSpace& pfes_scalar, - int vdim, const Vector& y_vdim_true, Vector& y) const; + void LumpedMassInverse(Vector& ML_inv) const override; + + const ParFiniteElementSpace& pfes_ho; + const ParFiniteElementSpace& pfes_lor; }; + #endif /** Mass-conservative prolongation operator going in the opposite direction diff --git a/miniapps/tools/CMakeLists.txt b/miniapps/tools/CMakeLists.txt index 36b612a702..cf6ba6966c 100644 --- a/miniapps/tools/CMakeLists.txt +++ b/miniapps/tools/CMakeLists.txt @@ -29,8 +29,8 @@ add_mfem_miniapp(lor-transfer MAIN lor-transfer.cpp LIBRARIES mfem) if (MFEM_USE_MPI) - add_mfem_miniapp(lor-transferp - MAIN lor-transferp.cpp LIBRARIES mfem) + add_mfem_miniapp(plor-transfer + MAIN plor-transfer.cpp LIBRARIES mfem) endif() add_mfem_miniapp(check-tmop-metric diff --git a/miniapps/tools/lor-transfer.cpp b/miniapps/tools/lor-transfer.cpp index 97695a07e8..38c2eaec36 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); 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/lor-transferp.cpp b/miniapps/tools/plor-transfer.cpp similarity index 87% rename from miniapps/tools/lor-transferp.cpp rename to miniapps/tools/plor-transfer.cpp index cabe13ddf1..21fb9ede58 100644 --- a/miniapps/tools/lor-transferp.cpp +++ b/miniapps/tools/plor-transfer.cpp @@ -29,17 +29,17 @@ // particular finite element spaces. For example they satisfy PR=I, plus mass // conservation in both directions for L2 fields. // -// Compile with: make lor-transferp +// Compile with: make plor-transfer // -// Sample runs: lor-transferp -// lor-transferp -h1 -// lor-transferp -t -// lor-transferp -m ../../data/star-q2.mesh -lref 5 -p 4 -// lor-transferp -m ../../data/star-mixed.mesh -lref 3 -p 2 -// lor-transferp -lref 4 -o 4 -lo 0 -p 1 -// lor-transferp -lref 5 -o 4 -lo 0 -p 1 -// lor-transferp -lref 5 -o 4 -lo 3 -p 2 -// lor-transferp -lref 5 -o 4 -lo 0 -p 3 +// 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 @@ -99,19 +99,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()) - { - if (Mpi::Root()) - { - args.PrintUsage(cout); - } - return 1; - } - if (Mpi::Root()) - { - args.PrintOptions(cout); - } + args.ParseCheck(); // Read the mesh from the given mesh file. Mesh serial_mesh(mesh_file, 1, 1); @@ -193,6 +181,11 @@ int main(int argc, char *argv[]) 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()) { @@ -207,9 +200,7 @@ int main(int argc, char *argv[]) rho_prev -= rho; Vector rho_prev_true(fespace.GetTrueVSize()); rho_prev.GetTrueDofs(rho_prev_true); - double l_inf_local = rho_prev_true.Normlinf(); - double l_inf; - MPI_Allreduce(&l_inf_local, &l_inf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + double l_inf = global_max(rho_prev_true); if (Mpi::Root()) { cout.precision(12); @@ -222,6 +213,11 @@ int main(int argc, char *argv[]) ones = 1.0; ones_lor = 1.0; 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()) { const Operator &P = gt->BackwardOperator(); @@ -234,13 +230,8 @@ int main(int argc, char *argv[]) Vector M_rho_lor_true(M_rho_lor.ParFESpace()->GetTrueVSize()); M_rho_lor.ParFESpace()->GetRestrictionOperator()->Mult(M_rho_lor, M_rho_lor_true); - double local_ho_mass = M_rho_true.Sum(); - double local_lor_mass = M_rho_lor_true.Sum(); - double ho_mass; - double lor_mass; - MPI_Allreduce(&local_ho_mass, &ho_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&local_lor_mass, &lor_mass, 1, MPI_DOUBLE, MPI_SUM, - MPI_COMM_WORLD); + double ho_mass = global_sum(M_rho_true); + double lor_mass = global_sum(M_rho_lor_true); if (Mpi::Root()) { cout << "HO -> LOR dual field: " << fabs(ho_mass - lor_mass) << endl << endl; @@ -273,9 +264,7 @@ int main(int argc, char *argv[]) rho_lor_prev -= rho_lor; Vector rho_lor_prev_true(fespace_lor.GetTrueVSize()); rho_lor_prev.GetTrueDofs(rho_lor_prev_true); - double l_inf_local = rho_lor_prev_true.Normlinf(); - double l_inf; - MPI_Allreduce(&l_inf_local, &l_inf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + double l_inf = global_max(rho_lor_prev_true); if (Mpi::Root()) { cout.precision(12); @@ -295,21 +284,14 @@ int main(int argc, char *argv[]) R.MultTranspose(M_rho_lor, M_rho); Vector M_rho_true(M_rho.ParFESpace()->GetTrueVSize()); M_rho.ParFESpace()->GetRestrictionOperator()->Mult(M_rho, M_rho_true); - double local_ho_mass = M_rho_true.Sum(); - double local_lor_mass = M_rho_lor_true.Sum(); - double ho_mass; - double lor_mass; - MPI_Allreduce(&local_ho_mass, &ho_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&local_lor_mass, &lor_mass, 1, MPI_DOUBLE, MPI_SUM, - MPI_COMM_WORLD); + double ho_mass = global_sum(M_rho_true); + double lor_mass = global_sum(M_rho_lor_true); if (Mpi::Root()) { cout << "LOR -> HO dual field: " << fabs(ho_mass - lor_mass) << '\n'; } } - Mpi::Finalize(); - delete fec; delete fec_lor; delete M_ho_tdof; From f12137346717700ecccca897b49f01bdca892aa5 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Sat, 29 Jul 2023 15:30:51 -0700 Subject: [PATCH 13/34] bugfixes --- fem/transfer.cpp | 28 +++++++++++++++++++--------- fem/transfer.hpp | 2 +- 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 677bfbc120..89af1baa62 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -585,7 +585,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( SetTDofsByVDim(fes_lor, d, Y_dim, Y); } - SetTDofs(fes_lor, Y, y); + SetTDofs(fes_lor, std::move(Y), y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( @@ -606,7 +606,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( SetTDofsByVDim(fes_ho, d, Y_dim, Y); } - SetTDofs(fes_ho, Y, y); + SetTDofs(fes_ho, std::move(Y), y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( @@ -631,7 +631,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( SetTDofsByVDim(fes_ho, d, Y_dim, Y); } - SetTDofs(fes_ho, Y, y); + SetTDofs(fes_ho, std::move(Y), y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( @@ -653,10 +653,10 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( Xbar = 0.0; pcg.Mult(X_dim, Xbar); M_LH->Mult(Xbar, Y_dim); - SetTDofsByVDim(fes_lor, d, Y_dim, y); + SetTDofsByVDim(fes_lor, d, Y_dim, Y); } - SetTDofs(fes_lor, Y, y); + SetTDofs(fes_lor, std::move(Y), y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_) @@ -747,9 +747,19 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() // 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) + 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( - r_and_mlh.first->GetI(), r_and_mlh.first->GetJ(), NULL, - r_and_mlh.first->Height(), r_and_mlh.first->Width(), false, true, true)); + 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; @@ -817,7 +827,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofs( } void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofs( - const FiniteElementSpace& fes, const Vector& X, Vector& x) const + const FiniteElementSpace& fes, Vector&& X, Vector& x) const { const Operator* P = fes.GetProlongationMatrix(); if (P) @@ -826,7 +836,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofs( } else { - x.NewDataAndSize(X.GetData(), X.Size()); + x = std::move(X); } } diff --git a/fem/transfer.hpp b/fem/transfer.hpp index f6d7ef3bb2..7e19b18e99 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -302,7 +302,7 @@ protected: virtual void GetTDofs(const FiniteElementSpace& fes, const Vector& x, Vector& X) const; virtual void SetTDofs(const FiniteElementSpace& fes, - const Vector& X, Vector& x) const; + Vector&& X, Vector& x) const; virtual void GetTDofsByVDim(const FiniteElementSpace& fes, int vdim, const Vector& X, Vector& X_vdim) const; virtual void SetTDofsByVDim(const FiniteElementSpace& fes, From fdd1c6c8b4c8c9472c845d4db46b49127a090bc0 Mon Sep 17 00:00:00 2001 From: EB Chin Date: Tue, 1 Aug 2023 13:20:31 -0700 Subject: [PATCH 14/34] simplify tdof method --- fem/transfer.cpp | 102 ++++++++++++++++++++++------------------------- fem/transfer.hpp | 15 +++---- 2 files changed, 54 insertions(+), 63 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 89af1baa62..c2363e5f59 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -543,7 +543,8 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::ProlongateTranspose( } L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( - const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_) + const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_, + bool build_operators) : L2Projection(fes_ho_, fes_lor_), #ifdef MFEM_USE_MPI pcg { dynamic_cast(&fes_ho) ? @@ -552,9 +553,12 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( pcg() #endif { - std::tie(R, M_LH) = ComputeSparseRAndM_LH(); - RTxM_LH = std::unique_ptr(TransposeMult(static_cast(*R), - static_cast(*M_LH))); + if (build_operators) + { + std::tie(R, M_LH) = ComputeSparseRAndM_LH(); + RTxM_LH = std::unique_ptr(TransposeMult(static_cast(*R), + static_cast(*M_LH))); + } // Basic PCG solver setup pcg.SetPrintLevel(0); @@ -562,9 +566,12 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( // initial values for relative and absolute tolerance SetRelTol(1e-13); SetAbsTol(1e-13); - precon = std::unique_ptr(new DSmoother(static_cast(*RTxM_LH))); - pcg.SetPreconditioner(*precon); - pcg.SetOperator(*RTxM_LH); + if (build_operators) + { + precon = std::unique_ptr(new DSmoother(static_cast(*RTxM_LH))); + pcg.SetPreconditioner(*precon); + pcg.SetOperator(*RTxM_LH); + } } void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( @@ -575,14 +582,18 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( 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) { - GetTDofsByVDim(fes_ho, d, X, X_dim); + TDofsListByVDim(fes_ho, d, vdofs_list); + X.GetSubVector(vdofs_list, X_dim); R->Mult(X_dim, Y_dim); - SetTDofsByVDim(fes_lor, d, Y_dim, Y); + TDofsListByVDim(fes_lor, d, vdofs_list); + Y.SetSubVector(vdofs_list, Y_dim); } SetTDofs(fes_lor, std::move(Y), y); @@ -597,13 +608,17 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( Vector Y_dim(R->Width()); Vector Y(fes_ho.GetTrueVSize()); + Array vdofs_list; + GetTDofs(fes_lor, x, X); for (int d = 0; d < fes_ho.GetVDim(); ++d) { - GetTDofsByVDim(fes_lor, d, X, X_dim); + TDofsListByVDim(fes_lor, d, vdofs_list); + X.GetSubVector(vdofs_list, X_dim); R->MultTranspose(X_dim, Y_dim); - SetTDofsByVDim(fes_ho, d, Y_dim, Y); + TDofsListByVDim(fes_ho, d, vdofs_list); + Y.SetSubVector(vdofs_list, Y_dim); } SetTDofs(fes_ho, std::move(Y), y); @@ -619,16 +634,20 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( 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) { - GetTDofsByVDim(fes_lor, d, X, X_dim); + 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); - SetTDofsByVDim(fes_ho, d, Y_dim, Y); + TDofsListByVDim(fes_ho, d, vdofs_list); + Y.SetSubVector(vdofs_list, Y_dim); } SetTDofs(fes_ho, std::move(Y), y); @@ -643,17 +662,21 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( Vector Y_dim(M_LH->Height()); Vector Y(fes_lor.GetTrueVSize()); + + Array vdofs_list; GetTDofs(fes_ho, x, X); for (int d = 0; d < fes_ho.GetVDim(); ++d) { - GetTDofsByVDim(fes_ho, d, X, X_dim); + 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); - SetTDofsByVDim(fes_lor, d, Y_dim, Y); + TDofsListByVDim(fes_lor, d, vdofs_list); + Y.SetSubVector(vdofs_list, Y_dim); } SetTDofs(fes_lor, std::move(Y), y); @@ -840,22 +863,11 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofs( } } -void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofsByVDim( - const FiniteElementSpace& fes, int vdim, const Vector& X, Vector& X_vdim) const +void L2ProjectionGridTransfer::L2ProjectionH1Space::TDofsListByVDim( + const FiniteElementSpace& fes, int vdim, Array& vdofs_list) const { - Array x_vdofs_list(fes.GetNDofs()); - - fes.GetVDofs(vdim, x_vdofs_list); - X.GetSubVector(x_vdofs_list, X_vdim); -} - -void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofsByVDim( - const FiniteElementSpace& fes, int vdim, const Vector& X_vdim, Vector& X) const -{ - Array x_vdofs_list(fes.GetNDofs()); - - fes.GetVDofs(vdim, x_vdofs_list); - X.SetSubVector(x_vdofs_list, X_vdim); + vdofs_list.SetSize(fes.GetNDofs()); + fes.GetVDofs(vdim, vdofs_list); } void L2ProjectionGridTransfer::L2ProjectionH1Space::LumpedMassInverse( @@ -960,10 +972,12 @@ std::unique_ptr L2ProjectionGridTransfer::L2ProjectionH1Space::All L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( const ParFiniteElementSpace& pfes_ho_, const ParFiniteElementSpace& pfes_lor_) - : L2ProjectionH1Space(pfes_ho_, pfes_lor_), + : L2ProjectionH1Space(pfes_ho_, pfes_lor_, false), pfes_ho(pfes_ho_), pfes_lor(pfes_lor_) { + 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(), @@ -998,37 +1012,17 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( pcg.SetOperator(*RTxM_LH); } -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::GetTDofsByVDim( - const FiniteElementSpace& fes, int vdim, const Vector& X, Vector& X_vdim) const +void L2ProjectionGridTransfer::ParL2ProjectionH1Space::TDofsListByVDim( + const FiniteElementSpace& fes, int vdim, Array& vdofs_list) const { Array x_vdofs_list(fes.GetNDofs()); Array x_vdofs_marker(fes.GetVSize()); Array X_vdofs_marker(fes.GetTrueVSize()); - Array X_vdofs_list; fes.GetVDofs(vdim, x_vdofs_list); FiniteElementSpace::ListToMarker(x_vdofs_list, fes.GetVSize(), x_vdofs_marker); fes.GetRestrictionMatrix()->BooleanMult(x_vdofs_marker, X_vdofs_marker); - FiniteElementSpace::MarkerToList(X_vdofs_marker, X_vdofs_list); - X.GetSubVector(X_vdofs_list, X_vdim); -} - -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::SetTDofsByVDim( - const FiniteElementSpace& fes, int vdim, const Vector& X_vdim, Vector& X) const -{ - const ParFiniteElementSpace* pfes = dynamic_cast(&fes); - MFEM_ASSERT(pfes != nullptr, "fes must be a parallel space"); - - Array x_vdofs_list(fes.GetNDofs()); - Array x_vdofs_marker(fes.GetVSize()); - Array X_vdofs_marker(fes.GetTrueVSize()); - Array X_vdofs_list; - - fes.GetVDofs(vdim, x_vdofs_list); - FiniteElementSpace::ListToMarker(x_vdofs_list, fes.GetVSize(), x_vdofs_marker); - pfes->Dof_TrueDof_Matrix()->BooleanMult(1.0, x_vdofs_marker, 0.0, X_vdofs_marker); - FiniteElementSpace::MarkerToList(X_vdofs_marker, X_vdofs_list); - X.SetSubVector(X_vdofs_list, X_vdim); + FiniteElementSpace::MarkerToList(X_vdofs_marker, vdofs_list); } void L2ProjectionGridTransfer::ParL2ProjectionH1Space::LumpedMassInverse( diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 7e19b18e99..35157a1a19 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -260,7 +260,8 @@ protected: { public: L2ProjectionH1Space(const FiniteElementSpace& fes_ho_, - const FiniteElementSpace& fes_lor_); + const FiniteElementSpace& fes_lor_, + bool build_operators = true); virtual ~L2ProjectionH1Space() = default; /// Maps x, primal field coefficients defined on a coarse mesh /// with a higher order H1 finite element space, to y, primal @@ -303,10 +304,8 @@ protected: const Vector& x, Vector& X) const; virtual void SetTDofs(const FiniteElementSpace& fes, Vector&& X, Vector& x) const; - virtual void GetTDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& X, Vector& X_vdim) const; - virtual void SetTDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& X_vdim, Vector& X) const; + virtual void TDofsListByVDim(const FiniteElementSpace& fes, + int vdim, Array& vdofs_list) const; /// Returns the inverse of an on-rank lumped mass matrix virtual void LumpedMassInverse(Vector& ML_inv) const; @@ -339,10 +338,8 @@ protected: const ParFiniteElementSpace& pfes_lor_); virtual ~ParL2ProjectionH1Space() = default; private: - void GetTDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& x, Vector& x_vdim) const override; - void SetTDofsByVDim(const FiniteElementSpace& fes, - int vdim, const Vector& x_vdim, Vector& x) const override; + void TDofsListByVDim(const FiniteElementSpace& fes, + int vdim, Array& vdofs_list) const override; /// Computes inverse of a lumped mass matrix (stored as a vector) void LumpedMassInverse(Vector& ML_inv) const override; From c0dab3037527fd551100c7bde0b2040eb68d2203 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Tue, 1 Aug 2023 22:00:26 -0700 Subject: [PATCH 15/34] rename variables --- miniapps/tools/plor-transfer.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/miniapps/tools/plor-transfer.cpp b/miniapps/tools/plor-transfer.cpp index 21fb9ede58..8cd5a578af 100644 --- a/miniapps/tools/plor-transfer.cpp +++ b/miniapps/tools/plor-transfer.cpp @@ -230,11 +230,11 @@ int main(int argc, char *argv[]) Vector M_rho_lor_true(M_rho_lor.ParFESpace()->GetTrueVSize()); M_rho_lor.ParFESpace()->GetRestrictionOperator()->Mult(M_rho_lor, M_rho_lor_true); - double ho_mass = global_sum(M_rho_true); - double lor_mass = global_sum(M_rho_lor_true); + double ho_dual_mass = global_sum(M_rho_true); + double lor_dual_mass = global_sum(M_rho_lor_true); if (Mpi::Root()) { - cout << "HO -> LOR dual field: " << fabs(ho_mass - lor_mass) << endl << endl; + cout << "HO -> LOR dual field: " << fabs(ho_dual_mass - lor_dual_mass) << endl << endl; } } @@ -284,11 +284,11 @@ int main(int argc, char *argv[]) R.MultTranspose(M_rho_lor, M_rho); Vector M_rho_true(M_rho.ParFESpace()->GetTrueVSize()); M_rho.ParFESpace()->GetRestrictionOperator()->Mult(M_rho, M_rho_true); - double ho_mass = global_sum(M_rho_true); - double lor_mass = global_sum(M_rho_lor_true); + double ho_dual_mass = global_sum(M_rho_true); + double lor_dual_mass = global_sum(M_rho_lor_true); if (Mpi::Root()) { - cout << "LOR -> HO dual field: " << fabs(ho_mass - lor_mass) << '\n'; + cout << "LOR -> HO dual field: " << fabs(ho_dual_mass - lor_dual_mass) << '\n'; } } From 7d2e402a08678afd6be3eb5af9b4d31ca2e1699f Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Tue, 1 Aug 2023 22:06:07 -0700 Subject: [PATCH 16/34] fix style --- fem/transfer.cpp | 46 +++++++++++++++++--------------- fem/transfer.hpp | 2 +- miniapps/tools/plor-transfer.cpp | 9 ++++--- 3 files changed, 32 insertions(+), 25 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index c2363e5f59..180fd55542 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -547,8 +547,9 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( bool build_operators) : L2Projection(fes_ho_, fes_lor_), #ifdef MFEM_USE_MPI - pcg { dynamic_cast(&fes_ho) ? - CGSolver(MPI_COMM_WORLD) : CGSolver() } + pcg { dynamic_cast(&fes_ho) ? + CGSolver(MPI_COMM_WORLD) : CGSolver() + } #else pcg() #endif @@ -556,10 +557,11 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( if (build_operators) { std::tie(R, M_LH) = ComputeSparseRAndM_LH(); - RTxM_LH = std::unique_ptr(TransposeMult(static_cast(*R), - static_cast(*M_LH))); + RTxM_LH = std::unique_ptr(TransposeMult(static_cast + (*R), + static_cast(*M_LH))); } - + // Basic PCG solver setup pcg.SetPrintLevel(0); pcg.SetMaxIter(1000); @@ -568,7 +570,8 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( SetAbsTol(1e-13); if (build_operators) { - precon = std::unique_ptr(new DSmoother(static_cast(*RTxM_LH))); + precon = std::unique_ptr(new DSmoother(static_cast + (*RTxM_LH))); pcg.SetPreconditioner(*precon); pcg.SetOperator(*RTxM_LH); } @@ -579,12 +582,12 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( { 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) @@ -664,7 +667,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( Vector Y(fes_lor.GetTrueVSize()); Array vdofs_list; - + GetTDofs(fes_ho, x, X); for (int d = 0; d < fes_ho.GetVDim(); ++d) @@ -693,10 +696,10 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_) } std::pair, std::unique_ptr> -L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() + L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() { std::pair, std::unique_ptr> - r_and_mlh; + r_and_mlh; Mesh* mesh_ho = fes_ho.GetMesh(); Mesh* mesh_lor = fes_lor.GetMesh(); @@ -781,8 +784,8 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() 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)); + 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; @@ -879,7 +882,8 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::LumpedMassInverse( } } -std::unique_ptr 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(); @@ -960,7 +964,7 @@ std::unique_ptr L2ProjectionGridTransfer::L2ProjectionH1Space::All double* data = Memory(dof_dofI[ndof_lor]); std::unique_ptr R_local(new SparseMatrix( - dof_dofI, dof_dofJ, data, ndof_lor, ndof_ho, true, true, true)); + dof_dofI, dof_dofJ, data, ndof_lor, ndof_ho, true, true, true)); (*R_local) = 0.0; dof_lor_dof_ho.LoseData(); @@ -993,20 +997,20 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( pfes_lor_scalar.GlobalVSize(), pfes_ho_scalar.GlobalVSize(), pfes_lor_scalar.GetDofOffsets(), - pfes_ho_scalar.GetDofOffsets(), + pfes_ho_scalar.GetDofOffsets(), static_cast(M_LH.get())); R = std::unique_ptr(RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), - &R_local, pfes_ho_scalar.Dof_TrueDof_Matrix())); + &R_local, pfes_ho_scalar.Dof_TrueDof_Matrix())); M_LH = std::unique_ptr(RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), - &M_LH_local, pfes_ho_scalar.Dof_TrueDof_Matrix())); - std::unique_ptr R_T = + &M_LH_local, pfes_ho_scalar.Dof_TrueDof_Matrix())); + std::unique_ptr R_T = std::unique_ptr(static_cast(*R).Transpose()); RTxM_LH = std::unique_ptr(ParMult( - R_T.get(), static_cast(M_LH.get()), true)); + R_T.get(), static_cast(M_LH.get()), true)); precon = std::unique_ptr( - new HypreBoomerAMG(static_cast(*RTxM_LH))); + new HypreBoomerAMG(static_cast(*RTxM_LH))); static_cast(*precon).SetPrintLevel(0); pcg.SetPreconditioner(*precon); pcg.SetOperator(*RTxM_LH); diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 35157a1a19..f558e6d81d 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -299,7 +299,7 @@ protected: protected: /// Computes on-rank R and M_LH matrices std::pair, std::unique_ptr> - ComputeSparseRAndM_LH(); + ComputeSparseRAndM_LH(); virtual void GetTDofs(const FiniteElementSpace& fes, const Vector& x, Vector& X) const; virtual void SetTDofs(const FiniteElementSpace& fes, diff --git a/miniapps/tools/plor-transfer.cpp b/miniapps/tools/plor-transfer.cpp index 8cd5a578af..a7dd4d4597 100644 --- a/miniapps/tools/plor-transfer.cpp +++ b/miniapps/tools/plor-transfer.cpp @@ -181,7 +181,8 @@ int main(int argc, char *argv[]) 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){ + 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; @@ -213,7 +214,8 @@ int main(int argc, char *argv[]) ones = 1.0; ones_lor = 1.0; ParLinearForm M_rho(&fespace), M_rho_lor(&fespace_lor); - auto global_sum = [](const Vector& v){ + 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; @@ -234,7 +236,8 @@ int main(int argc, char *argv[]) double lor_dual_mass = global_sum(M_rho_lor_true); if (Mpi::Root()) { - cout << "HO -> LOR dual field: " << fabs(ho_dual_mass - lor_dual_mass) << endl << endl; + cout << "HO -> LOR dual field: " << fabs(ho_dual_mass - lor_dual_mass) << endl + << endl; } } From 100de86cfa0c28f9824acebea5ade0e875c4b12d Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Tue, 1 Aug 2023 22:14:22 -0700 Subject: [PATCH 17/34] add new miniapp and new feature --- CHANGELOG | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/CHANGELOG b/CHANGELOG index db8260bd68..666e488de1 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -39,6 +39,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 new methods in the Mesh class to set and get attributes on NURBS patches @@ -59,6 +63,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 From e42147f3a7a3f9eb017a9ff7324ea519cf05ad86 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Mon, 7 Aug 2023 11:14:40 -0700 Subject: [PATCH 18/34] Simplify ParL2ProjectionH1Space --- fem/transfer.cpp | 105 ++++++++++++++++++++++------------------------- fem/transfer.hpp | 29 ++++--------- 2 files changed, 59 insertions(+), 75 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 180fd55542..7d31206ea1 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -869,17 +869,46 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofs( void L2ProjectionGridTransfer::L2ProjectionH1Space::TDofsListByVDim( const FiniteElementSpace& fes, int vdim, Array& vdofs_list) const { - vdofs_list.SetSize(fes.GetNDofs()); - fes.GetVDofs(vdim, vdofs_list); + 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 + { + fes.GetVDofs(vdim, vdofs_list); + } } void L2ProjectionGridTransfer::L2ProjectionH1Space::LumpedMassInverse( Vector& ML_inv) const { - for (int i = 0; i < ML_inv.Size(); ++i) + 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[i] = 1.0 / ML_inv[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); } std::unique_ptr @@ -975,10 +1004,8 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() #ifdef MFEM_USE_MPI L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( - const ParFiniteElementSpace& pfes_ho_, const ParFiniteElementSpace& pfes_lor_) - : L2ProjectionH1Space(pfes_ho_, pfes_lor_, false), - pfes_ho(pfes_ho_), - pfes_lor(pfes_lor_) + const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor) + : L2ProjectionH1Space(pfes_ho, pfes_lor, false) { std::tie(R, M_LH) = ComputeSparseRAndM_LH(); @@ -1000,58 +1027,26 @@ L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( pfes_ho_scalar.GetDofOffsets(), static_cast(M_LH.get())); - R = std::unique_ptr(RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), - &R_local, pfes_ho_scalar.Dof_TrueDof_Matrix())); - M_LH = std::unique_ptr(RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(), - &M_LH_local, pfes_ho_scalar.Dof_TrueDof_Matrix())); - std::unique_ptr R_T = - std::unique_ptr(static_cast(*R).Transpose()); - RTxM_LH = std::unique_ptr(ParMult( - R_T.get(), static_cast(M_LH.get()), true)); + 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); - precon = std::unique_ptr( - new HypreBoomerAMG(static_cast(*RTxM_LH))); - static_cast(*precon).SetPrintLevel(0); pcg.SetPreconditioner(*precon); pcg.SetOperator(*RTxM_LH); } -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::TDofsListByVDim( - const FiniteElementSpace& fes, int vdim, Array& vdofs_list) const -{ - 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); - fes.GetRestrictionMatrix()->BooleanMult(x_vdofs_marker, X_vdofs_marker); - FiniteElementSpace::MarkerToList(X_vdofs_marker, vdofs_list); -} - -void L2ProjectionGridTransfer::ParL2ProjectionH1Space::LumpedMassInverse( - Vector& ML_inv) const -{ - Vector ML_inv_full(pfes_lor.GetVSize()); - // set ML_inv on dofs for vdim = 0 - Array vdofs_list(pfes_lor.GetNDofs()); - pfes_lor.GetVDofs(0, vdofs_list); - ML_inv_full.SetSubVector(vdofs_list, ML_inv); - - Vector ML_inv_true(pfes_lor.GetTrueVSize()); - const Operator& P = *pfes_lor.GetProlongationMatrix(); - P.MultTranspose(ML_inv_full, ML_inv_true); - - for (int i = 0; i < ML_inv_true.Size(); ++i) - { - ML_inv_true[i] = 1.0 / ML_inv_true[i]; - } - - P.Mult(ML_inv_true, ML_inv_full); - - ML_inv_full.GetSubVector(vdofs_list, ML_inv); -} - #endif L2ProjectionGridTransfer::~L2ProjectionGridTransfer() diff --git a/fem/transfer.hpp b/fem/transfer.hpp index f558e6d81d..35a7f84c4a 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -298,16 +298,14 @@ protected: virtual void SetAbsTol(double p_atol_); protected: /// Computes on-rank R and M_LH matrices - std::pair, std::unique_ptr> - ComputeSparseRAndM_LH(); - virtual void GetTDofs(const FiniteElementSpace& fes, - const Vector& x, Vector& X) const; - virtual void SetTDofs(const FiniteElementSpace& fes, - Vector&& X, Vector& x) const; - virtual void TDofsListByVDim(const FiniteElementSpace& fes, - int vdim, Array& vdofs_list) const; + std::pair, + std::unique_ptr> ComputeSparseRAndM_LH(); + void GetTDofs(const FiniteElementSpace& fes, const Vector& x, Vector& X) const; + void SetTDofs(const FiniteElementSpace& fes, Vector&& X, Vector& x) const; + void TDofsListByVDim(const FiniteElementSpace& fes, + int vdim, Array& vdofs_list) const; /// Returns the inverse of an on-rank lumped mass matrix - virtual void LumpedMassInverse(Vector& ML_inv) const; + void LumpedMassInverse(Vector& ML_inv) const; CGSolver pcg; std::unique_ptr precon; @@ -334,17 +332,8 @@ protected: class ParL2ProjectionH1Space : public L2ProjectionH1Space { public: - ParL2ProjectionH1Space(const ParFiniteElementSpace& pfes_ho_, - const ParFiniteElementSpace& pfes_lor_); - virtual ~ParL2ProjectionH1Space() = default; - private: - void TDofsListByVDim(const FiniteElementSpace& fes, - int vdim, Array& vdofs_list) const override; - /// Computes inverse of a lumped mass matrix (stored as a vector) - void LumpedMassInverse(Vector& ML_inv) const override; - - const ParFiniteElementSpace& pfes_ho; - const ParFiniteElementSpace& pfes_lor; + ParL2ProjectionH1Space(const ParFiniteElementSpace& pfes_ho, + const ParFiniteElementSpace& pfes_lor); }; #endif From 39714f039a9996a3f6018a9f1e5ab96a29b78c9d Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Mon, 7 Aug 2023 11:34:16 -0700 Subject: [PATCH 19/34] Merge ParL2ProjectionH1Space with L2ProjectionH1Space --- fem/transfer.cpp | 180 +++++++++++++++++++++++++++-------------------- fem/transfer.hpp | 26 +++---- 2 files changed, 111 insertions(+), 95 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 7d31206ea1..8331c23bd6 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -543,38 +543,109 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::ProlongateTranspose( } L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space( - const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_, - bool build_operators) - : L2Projection(fes_ho_, fes_lor_), -#ifdef MFEM_USE_MPI - pcg { dynamic_cast(&fes_ho) ? - CGSolver(MPI_COMM_WORLD) : CGSolver() - } -#else - pcg() -#endif + const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_) + : L2Projection(fes_ho_, fes_lor_) { - if (build_operators) + 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) { - std::tie(R, M_LH) = ComputeSparseRAndM_LH(); - RTxM_LH = std::unique_ptr(TransposeMult(static_cast - (*R), - static_cast(*M_LH))); + 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 - SetRelTol(1e-13); - SetAbsTol(1e-13); - if (build_operators) - { - precon = std::unique_ptr(new DSmoother(static_cast - (*RTxM_LH))); - pcg.SetPreconditioner(*precon); - pcg.SetOperator(*RTxM_LH); - } + pcg.SetRelTol(1e-13); + pcg.SetAbsTol(1e-13); + pcg.SetPreconditioner(*precon); + pcg.SetOperator(*RTxM_LH); } void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( @@ -695,11 +766,13 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_) pcg.SetAbsTol(p_atol_); } -std::pair, std::unique_ptr> - L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() +std::pair< +std::unique_ptr, +std::unique_ptr> + L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH() { - std::pair, std::unique_ptr> - r_and_mlh; + std::pair, + std::unique_ptr> r_and_mlh; Mesh* mesh_ho = fes_ho.GetMesh(); Mesh* mesh_lor = fes_lor.GetMesh(); @@ -882,6 +955,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::TDofsListByVDim( } else { + vdofs_list.SetSize(fes.GetNDofs()); fes.GetVDofs(vdim, vdofs_list); } } @@ -1001,54 +1075,6 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR() return R_local; } -#ifdef MFEM_USE_MPI - -L2ProjectionGridTransfer::ParL2ProjectionH1Space::ParL2ProjectionH1Space( - const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor) - : L2ProjectionH1Space(pfes_ho, pfes_lor, false) -{ - 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); - - pcg.SetPreconditioner(*precon); - pcg.SetOperator(*RTxM_LH); -} - -#endif - L2ProjectionGridTransfer::~L2ProjectionGridTransfer() { delete F; @@ -1087,7 +1113,7 @@ void L2ProjectionGridTransfer::BuildF() static_cast(dom_fes); const mfem::ParFiniteElementSpace& ran_pfes = static_cast(ran_fes); - F = new ParL2ProjectionH1Space(dom_pfes, ran_pfes); + F = new L2ProjectionH1Space(dom_pfes, ran_pfes); #endif } } diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 35a7f84c4a..01969c246b 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -259,9 +259,12 @@ protected: class L2ProjectionH1Space : public L2Projection { public: - L2ProjectionH1Space(const FiniteElementSpace& fes_ho_, - const FiniteElementSpace& fes_lor_, - bool build_operators = true); + L2ProjectionH1Space(const FiniteElementSpace &fes_ho_, + const FiniteElementSpace &fes_lor_); +#ifdef MFEM_USE_MPI + L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_, + const ParFiniteElementSpace &pfes_lor_); +#endif virtual ~L2ProjectionH1Space() = default; /// Maps x, primal field coefficients defined on a coarse mesh /// with a higher order H1 finite element space, to y, primal @@ -297,6 +300,8 @@ protected: virtual void SetRelTol(double p_rtol_); virtual void SetAbsTol(double p_atol_); 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(); @@ -323,21 +328,6 @@ protected: std::unique_ptr AllocR(); }; - -#ifdef MFEM_USE_MPI - - /** Implements 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) in parallel. */ - class ParL2ProjectionH1Space : public L2ProjectionH1Space - { - public: - ParL2ProjectionH1Space(const ParFiniteElementSpace& pfes_ho, - const ParFiniteElementSpace& pfes_lor); - }; - -#endif - /** Mass-conservative prolongation operator going in the opposite direction as L2Projection. This operator is a left inverse to the L2Projection. */ class L2Prolongation : public Operator From cddbd24df31f9b6bc5860bee7d2e91668b8424a7 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Mon, 7 Aug 2023 13:26:24 -0700 Subject: [PATCH 20/34] Simplify mass computation in lor-transfer miniapps --- miniapps/tools/lor-transfer.cpp | 11 ++++------- miniapps/tools/plor-transfer.cpp | 18 ++++-------------- 2 files changed, 8 insertions(+), 21 deletions(-) diff --git a/miniapps/tools/lor-transfer.cpp b/miniapps/tools/lor-transfer.cpp index 38c2eaec36..7c2269a411 100644 --- a/miniapps/tools/lor-transfer.cpp +++ b/miniapps/tools/lor-transfer.cpp @@ -282,14 +282,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/plor-transfer.cpp b/miniapps/tools/plor-transfer.cpp index a7dd4d4597..6dc776855e 100644 --- a/miniapps/tools/plor-transfer.cpp +++ b/miniapps/tools/plor-transfer.cpp @@ -345,21 +345,11 @@ double compute_mass(ParFiniteElementSpace *L2, double massL2, VisItDataCollection &dc, string prefix) { ConstantCoefficient one(1.0); - ParBilinearForm ML2(L2); - ML2.AddDomainIntegrator(new MassIntegrator(one)); - ML2.Assemble(); - ML2.Finalize(); - HypreParMatrix* pML2 = ML2.ParallelAssemble(); + ParLinearForm lf(L2); + lf.AddDomainIntegrator(new DomainLFIntegrator(one)); + lf.Assemble(); - Vector rhoone(L2->GetTrueVSize()); - rhoone = 1.0; - - Vector Mdiag(L2->GetTrueVSize()); - pML2->Mult(rhoone, Mdiag); - delete pML2; - HypreParVector* rho = dc.GetParField("density")->GetTrueDofs(); - double newmass = InnerProduct(MPI_COMM_WORLD, *rho, Mdiag); - delete rho; + double newmass = lf(*dc.GetParField("density")); if (Mpi::Root()) { cout.precision(18); From 89f7d276efaa6e73298e63cd4a64b11965fb284c Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Mon, 7 Aug 2023 13:26:32 -0700 Subject: [PATCH 21/34] Remove unused variables --- miniapps/tools/plor-transfer.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/miniapps/tools/plor-transfer.cpp b/miniapps/tools/plor-transfer.cpp index 6dc776855e..c4a534945a 100644 --- a/miniapps/tools/plor-transfer.cpp +++ b/miniapps/tools/plor-transfer.cpp @@ -210,9 +210,6 @@ int main(int argc, char *argv[]) } // HO* to LOR* dual fields - ParGridFunction ones(&fespace), ones_lor(&fespace_lor); - ones = 1.0; - ones_lor = 1.0; ParLinearForm M_rho(&fespace), M_rho_lor(&fespace_lor); auto global_sum = [](const Vector& v) { From c4dc57ffd242bb2a9d3dc1b0f1b298ab2b96a0ae Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Mon, 7 Aug 2023 13:34:08 -0700 Subject: [PATCH 22/34] Always perform copy in L2ProjectionH1Space GetTDofs and SetTDofs --- fem/transfer.cpp | 14 +++++++------- fem/transfer.hpp | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 8331c23bd6..9ba569a3c8 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -670,7 +670,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( Y.SetSubVector(vdofs_list, Y_dim); } - SetTDofs(fes_lor, std::move(Y), y); + SetTDofs(fes_lor, Y, y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( @@ -695,7 +695,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( Y.SetSubVector(vdofs_list, Y_dim); } - SetTDofs(fes_ho, std::move(Y), y); + SetTDofs(fes_ho, Y, y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( @@ -724,7 +724,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( Y.SetSubVector(vdofs_list, Y_dim); } - SetTDofs(fes_ho, std::move(Y), y); + SetTDofs(fes_ho, Y, y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( @@ -753,7 +753,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( Y.SetSubVector(vdofs_list, Y_dim); } - SetTDofs(fes_lor, std::move(Y), y); + SetTDofs(fes_lor, Y, y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_) @@ -921,12 +921,12 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofs( } else { - X.NewDataAndSize(x.GetData(), x.Size()); + X = x; } } void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofs( - const FiniteElementSpace& fes, Vector&& X, Vector& x) const + const FiniteElementSpace& fes, const Vector &X, Vector& x) const { const Operator* P = fes.GetProlongationMatrix(); if (P) @@ -935,7 +935,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofs( } else { - x = std::move(X); + x = X; } } diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 01969c246b..6863b2afcc 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -306,7 +306,7 @@ protected: std::pair, std::unique_ptr> ComputeSparseRAndM_LH(); void GetTDofs(const FiniteElementSpace& fes, const Vector& x, Vector& X) const; - void SetTDofs(const FiniteElementSpace& fes, Vector&& X, Vector& x) const; + void SetTDofs(const FiniteElementSpace& fes, const Vector& X, Vector& x) const; void TDofsListByVDim(const FiniteElementSpace& fes, int vdim, Array& vdofs_list) const; /// Returns the inverse of an on-rank lumped mass matrix From aeacd50a405ccced205cc3b922d5f0a61cd69c9b Mon Sep 17 00:00:00 2001 From: EB Chin Date: Tue, 8 Aug 2023 12:40:43 -0700 Subject: [PATCH 23/34] fix dual field parallel transfers --- fem/transfer.cpp | 42 ++++++++++++++++++++++++++------ fem/transfer.hpp | 10 +++++--- miniapps/tools/plor-transfer.cpp | 12 ++++----- 3 files changed, 47 insertions(+), 17 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 9ba569a3c8..4ec66ea633 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -670,7 +670,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult( Y.SetSubVector(vdofs_list, Y_dim); } - SetTDofs(fes_lor, Y, y); + SetFromTDofs(fes_lor, Y, y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( @@ -684,7 +684,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( Array vdofs_list; - GetTDofs(fes_lor, x, X); + GetTDofsTranspose(fes_lor, x, X); for (int d = 0; d < fes_ho.GetVDim(); ++d) { @@ -695,7 +695,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose( Y.SetSubVector(vdofs_list, Y_dim); } - SetTDofs(fes_ho, Y, y); + SetFromTDofsTranspose(fes_ho, Y, y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( @@ -724,7 +724,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate( Y.SetSubVector(vdofs_list, Y_dim); } - SetTDofs(fes_ho, Y, y); + SetFromTDofs(fes_ho, Y, y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( @@ -739,7 +739,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( Array vdofs_list; - GetTDofs(fes_ho, x, X); + GetTDofsTranspose(fes_ho, x, X); for (int d = 0; d < fes_ho.GetVDim(); ++d) { @@ -753,7 +753,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose( Y.SetSubVector(vdofs_list, Y_dim); } - SetTDofs(fes_lor, Y, y); + SetFromTDofsTranspose(fes_lor, Y, y); } void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_) @@ -925,7 +925,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofs( } } -void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofs( +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofs( const FiniteElementSpace& fes, const Vector &X, Vector& x) const { const Operator* P = fes.GetProlongationMatrix(); @@ -939,6 +939,34 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetTDofs( } } +void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofsTranspose( + const FiniteElementSpace& fes, const Vector& x, Vector& X) const +{ + const Operator* P = fes.GetProlongationMatrix(); + if (P) + { + P->MultTranspose(x, X); + } + else + { + X = x; + } +} + +void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofsTranspose( + const FiniteElementSpace& fes, const Vector &X, Vector& x) const +{ + const Operator* Rtranspose = fes.GetRestrictionTransposeOperator(); + if (Rtranspose) + { + Rtranspose->Mult(X, x); + } + else + { + x = X; + } +} + void L2ProjectionGridTransfer::L2ProjectionH1Space::TDofsListByVDim( const FiniteElementSpace& fes, int vdim, Array& vdofs_list) const { diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 6863b2afcc..7e679ed548 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -253,9 +253,9 @@ protected: virtual void SetAbsTol(double p_atol_) {} }; - /** Abstract 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 { public: @@ -306,7 +306,9 @@ protected: std::pair, std::unique_ptr> ComputeSparseRAndM_LH(); void GetTDofs(const FiniteElementSpace& fes, const Vector& x, Vector& X) const; - void SetTDofs(const FiniteElementSpace& fes, const Vector& X, Vector& x) const; + void SetFromTDofs(const FiniteElementSpace& fes, const Vector& X, Vector& x) const; + void GetTDofsTranspose(const FiniteElementSpace& fes, const Vector& x, Vector& X) const; + void SetFromTDofsTranspose(const FiniteElementSpace& fes, const Vector& X, Vector& x) const; void TDofsListByVDim(const FiniteElementSpace& fes, int vdim, Array& vdofs_list) const; /// Returns the inverse of an on-rank lumped mass matrix diff --git a/miniapps/tools/plor-transfer.cpp b/miniapps/tools/plor-transfer.cpp index c4a534945a..b0c3d00d41 100644 --- a/miniapps/tools/plor-transfer.cpp +++ b/miniapps/tools/plor-transfer.cpp @@ -224,11 +224,11 @@ int main(int argc, char *argv[]) rho.GetTrueDofs(rho_true); Vector M_rho_true(M_rho.ParFESpace()->GetTrueVSize()); M_ho_tdof->Mult(rho_true, M_rho_true); - M_rho.ParFESpace()->GetProlongationMatrix()->Mult(M_rho_true, M_rho); + M_rho.ParFESpace()->GetRestrictionTransposeOperator()->Mult(M_rho_true, M_rho); P.MultTranspose(M_rho, M_rho_lor); Vector M_rho_lor_true(M_rho_lor.ParFESpace()->GetTrueVSize()); - M_rho_lor.ParFESpace()->GetRestrictionOperator()->Mult(M_rho_lor, - M_rho_lor_true); + M_rho_lor.ParFESpace()->GetProlongationMatrix()->MultTranspose(M_rho_lor, + M_rho_lor_true); double ho_dual_mass = global_sum(M_rho_true); double lor_dual_mass = global_sum(M_rho_lor_true); if (Mpi::Root()) @@ -279,11 +279,11 @@ int main(int argc, char *argv[]) rho_lor.GetTrueDofs(rho_lor_true); Vector M_rho_lor_true(M_rho_lor.ParFESpace()->GetTrueVSize()); M_lor_tdof->Mult(rho_lor_true, M_rho_lor_true); - M_rho_lor.ParFESpace()->GetProlongationMatrix()->Mult(M_rho_lor_true, - M_rho_lor); + M_rho_lor.ParFESpace()->GetRestrictionTransposeOperator()->Mult(M_rho_lor_true, + M_rho_lor); R.MultTranspose(M_rho_lor, M_rho); Vector M_rho_true(M_rho.ParFESpace()->GetTrueVSize()); - M_rho.ParFESpace()->GetRestrictionOperator()->Mult(M_rho, M_rho_true); + M_rho.ParFESpace()->GetProlongationMatrix()->MultTranspose(M_rho, M_rho_true); double ho_dual_mass = global_sum(M_rho_true); double lor_dual_mass = global_sum(M_rho_lor_true); if (Mpi::Root()) From 3e1a24c35c772d1db45dce46f78b77816a6cfc47 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Tue, 8 Aug 2023 14:49:26 -0700 Subject: [PATCH 24/34] add amr option to lor-transfer --- miniapps/tools/lor-transfer.cpp | 202 ++++++++++++++++++-------------- 1 file changed, 114 insertions(+), 88 deletions(-) diff --git a/miniapps/tools/lor-transfer.cpp b/miniapps/tools/lor-transfer.cpp index 7c2269a411..fdcd08b098 100644 --- a/miniapps/tools/lor-transfer.cpp +++ b/miniapps/tools/lor-transfer.cpp @@ -75,6 +75,7 @@ int main(int argc, char *argv[]) bool vis = true; bool useH1 = false; bool use_pointwise_transfer = false; + int amr_iterations = 0; OptionsParser args(argc, argv); args.AddOption(&mesh_file, "-m", "--mesh", @@ -95,6 +96,8 @@ 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.AddOption(&amr_iterations, "-a", "--amr-iterations", + "Number of times to apply adaptive mesh refinement (AMR)."); args.ParseCheck(); // Read the mesh from the given mesh file. @@ -137,104 +140,127 @@ int main(int argc, char *argv[]) VisItDataCollection LOR_dc("LOR", &mesh_lor); LOR_dc.RegisterField("density", &rho_lor); - BilinearForm M_ho(&fespace); - M_ho.AddDomainIntegrator(new MassIntegrator); - M_ho.Assemble(); - M_ho.Finalize(); - - BilinearForm M_lor(&fespace_lor); - M_lor.AddDomainIntegrator(new MassIntegrator); - M_lor.Assemble(); - M_lor.Finalize(); - // HO projections direction = "HO -> LOR @ HO"; FunctionCoefficient RHO(RHO_exact); - rho.ProjectCoefficient(RHO); - 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) + CoefficientRefiner refiner(RHO, order); + refiner.SetThreshold(1.0e-4); + + for (int it = 0; it <= amr_iterations; ++it) { - gt = new InterpolationGridTransfer(fespace, fespace_lor); - } - else - { - gt = new L2ProjectionGridTransfer(fespace, fespace_lor); - } - const Operator &R = gt->ForwardOperator(); + rho.ProjectCoefficient(RHO); + + double ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO "); + if (vis) { visualize(HO_dc, "HO", Wx, Wy); Wx += offx; } - // 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; } + GridTransfer *gt; + if (use_pointwise_transfer) + { + gt = new InterpolationGridTransfer(fespace, fespace_lor); + } + else + { + gt = new L2ProjectionGridTransfer(fespace, fespace_lor); + } + const Operator &R = gt->ForwardOperator(); - if (gt->SupportsBackwardsOperator()) - { - const Operator &P = gt->BackwardOperator(); - // LOR->HO prolongation - direction = "HO -> LOR @ HO"; - GridFunction 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; - cout.precision(12); - cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl; - } - - // 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; - } - - // LOR projections - direction = "LOR -> HO @ LOR"; - rho_lor.ProjectCoefficient(RHO); - GridFunction 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"; + // HO->LOR restriction + direction = "HO -> LOR @ 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); } + compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO) "); + if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy); Wx += offx; } - rho_lor_prev -= rho_lor; - cout.precision(12); - cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl; - } + if (gt->SupportsBackwardsOperator()) + { + const Operator &P = gt->BackwardOperator(); + // LOR->HO prolongation + direction = "HO -> LOR @ HO"; + GridFunction 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; } - // LOR* to HO* dual fields - if (!use_pointwise_transfer) - { - 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'; + rho_prev -= rho; + cout.precision(12); + cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl; + } + + // 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(); + BilinearForm M_ho(&fespace); + M_ho.AddDomainIntegrator(new MassIntegrator); + M_ho.Assemble(); + M_ho.Finalize(); + 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; + } + + // LOR projections + direction = "LOR -> HO @ LOR"; + rho_lor.ProjectCoefficient(RHO); + GridFunction 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; + cout.precision(12); + cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl; + } + + // LOR* to HO* dual fields + if (!use_pointwise_transfer) + { + BilinearForm M_lor(&fespace_lor); + M_lor.AddDomainIntegrator(new MassIntegrator); + M_lor.Assemble(); + M_lor.Finalize(); + 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'; + } + + if (it < amr_iterations) + { + // call the refiner to modify the mesh + cout << "\nAMR refinement: Original element count = " << mesh.GetNE(); + refiner.Apply(mesh); + cout << " Post-AMR element count = " << mesh.GetNE() << endl; + + fespace.Update(); + rho.Update(); + + // build refined AMR modified mesh + mesh_lor = Mesh::MakeRefined(mesh, lref, basis_lor); + + fespace_lor.Update(); + rho_lor.Update(); + } } delete fec; From 32432c2c19bdce08f73c2c7ab182c1743e08b3a9 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Wed, 9 Aug 2023 11:16:43 -0700 Subject: [PATCH 25/34] Revert "add amr option to lor-transfer" This reverts commit 3e1a24c35c772d1db45dce46f78b77816a6cfc47. --- miniapps/tools/lor-transfer.cpp | 194 ++++++++++++++------------------ 1 file changed, 84 insertions(+), 110 deletions(-) diff --git a/miniapps/tools/lor-transfer.cpp b/miniapps/tools/lor-transfer.cpp index fdcd08b098..7c2269a411 100644 --- a/miniapps/tools/lor-transfer.cpp +++ b/miniapps/tools/lor-transfer.cpp @@ -75,7 +75,6 @@ int main(int argc, char *argv[]) bool vis = true; bool useH1 = false; bool use_pointwise_transfer = false; - int amr_iterations = 0; OptionsParser args(argc, argv); args.AddOption(&mesh_file, "-m", "--mesh", @@ -96,8 +95,6 @@ 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.AddOption(&amr_iterations, "-a", "--amr-iterations", - "Number of times to apply adaptive mesh refinement (AMR)."); args.ParseCheck(); // Read the mesh from the given mesh file. @@ -140,127 +137,104 @@ int main(int argc, char *argv[]) VisItDataCollection LOR_dc("LOR", &mesh_lor); LOR_dc.RegisterField("density", &rho_lor); + BilinearForm M_ho(&fespace); + M_ho.AddDomainIntegrator(new MassIntegrator); + M_ho.Assemble(); + M_ho.Finalize(); + + BilinearForm M_lor(&fespace_lor); + M_lor.AddDomainIntegrator(new MassIntegrator); + M_lor.Assemble(); + M_lor.Finalize(); + // HO projections direction = "HO -> LOR @ HO"; FunctionCoefficient RHO(RHO_exact); + rho.ProjectCoefficient(RHO); + double ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO "); + if (vis) { visualize(HO_dc, "HO", Wx, Wy); Wx += offx; } - CoefficientRefiner refiner(RHO, order); - refiner.SetThreshold(1.0e-4); - - for (int it = 0; it <= amr_iterations; ++it) + GridTransfer *gt; + if (use_pointwise_transfer) { - rho.ProjectCoefficient(RHO); - - double ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO "); - if (vis) { visualize(HO_dc, "HO", Wx, Wy); Wx += offx; } + gt = new InterpolationGridTransfer(fespace, fespace_lor); + } + else + { + gt = new L2ProjectionGridTransfer(fespace, fespace_lor); + } + const Operator &R = gt->ForwardOperator(); - 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; } - // 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; } + if (gt->SupportsBackwardsOperator()) + { + const Operator &P = gt->BackwardOperator(); + // LOR->HO prolongation + direction = "HO -> LOR @ HO"; + GridFunction 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; } - if (gt->SupportsBackwardsOperator()) - { - const Operator &P = gt->BackwardOperator(); - // LOR->HO prolongation - direction = "HO -> LOR @ HO"; - GridFunction 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; + cout.precision(12); + cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl; + } - rho_prev -= rho; - cout.precision(12); - cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl; - } + // 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; + } - // 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(); - BilinearForm M_ho(&fespace); - M_ho.AddDomainIntegrator(new MassIntegrator); - M_ho.Assemble(); - M_ho.Finalize(); - 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; - } + // LOR projections + direction = "LOR -> HO @ LOR"; + rho_lor.ProjectCoefficient(RHO); + GridFunction 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; } - // LOR projections + 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"; - rho_lor.ProjectCoefficient(RHO); - GridFunction 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; } + 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); } - 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; } + rho_lor_prev -= rho_lor; + cout.precision(12); + cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl; + } - // 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; - cout.precision(12); - cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl; - } - - // LOR* to HO* dual fields - if (!use_pointwise_transfer) - { - BilinearForm M_lor(&fespace_lor); - M_lor.AddDomainIntegrator(new MassIntegrator); - M_lor.Assemble(); - M_lor.Finalize(); - 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'; - } - - if (it < amr_iterations) - { - // call the refiner to modify the mesh - cout << "\nAMR refinement: Original element count = " << mesh.GetNE(); - refiner.Apply(mesh); - cout << " Post-AMR element count = " << mesh.GetNE() << endl; - - fespace.Update(); - rho.Update(); - - // build refined AMR modified mesh - mesh_lor = Mesh::MakeRefined(mesh, lref, basis_lor); - - fespace_lor.Update(); - rho_lor.Update(); - } + // LOR* to HO* dual fields + if (!use_pointwise_transfer) + { + 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'; } delete fec; From 1cf0ffbe0387d3dd5bab16e260cd0407d580aabc Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Wed, 9 Aug 2023 11:17:55 -0700 Subject: [PATCH 26/34] Remove virtual destructor from L2ProjectionH1Space --- fem/transfer.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 7e679ed548..2f8792bd3c 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -265,7 +265,6 @@ protected: L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_, const ParFiniteElementSpace &pfes_lor_); #endif - virtual ~L2ProjectionH1Space() = default; /// 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 From 65921a4ddab4903d6426eb866d059dd07694d581 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Wed, 9 Aug 2023 11:20:01 -0700 Subject: [PATCH 27/34] make style --- fem/transfer.hpp | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 2f8792bd3c..87cb0020e9 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -305,11 +305,18 @@ protected: std::pair, std::unique_ptr> ComputeSparseRAndM_LH(); void GetTDofs(const FiniteElementSpace& fes, const Vector& x, Vector& X) const; - void SetFromTDofs(const FiniteElementSpace& fes, const Vector& X, Vector& x) const; - void GetTDofsTranspose(const FiniteElementSpace& fes, const Vector& x, Vector& X) const; - void SetFromTDofsTranspose(const FiniteElementSpace& fes, const Vector& X, Vector& x) const; + void SetFromTDofs(const FiniteElementSpace& fes, + const Vector& X, + Vector& x) const; + void GetTDofsTranspose(const FiniteElementSpace& fes, + const Vector& x, + Vector& X) const; + void SetFromTDofsTranspose(const FiniteElementSpace& fes, + const Vector& X, + Vector& x) const; void TDofsListByVDim(const FiniteElementSpace& fes, - int vdim, Array& vdofs_list) const; + int vdim, + Array& vdofs_list) const; /// Returns the inverse of an on-rank lumped mass matrix void LumpedMassInverse(Vector& ML_inv) const; From 1a10217ba481b63bfb6fad77c3ae2c82a257daf6 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Wed, 9 Aug 2023 11:20:21 -0700 Subject: [PATCH 28/34] Move L2ProjectionH1Space::AllocR from private to protected --- fem/transfer.hpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 87cb0020e9..2745e8b092 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -319,6 +319,10 @@ protected: Array& vdofs_list) const; /// Returns the inverse of an on-rank lumped mass matrix void LumpedMassInverse(Vector& ML_inv) const; + /// 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; @@ -329,11 +333,6 @@ protected: // Used to compute P = (RT*M_LH)^(-1) M_LH^T std::unique_ptr M_LH; std::unique_ptr RTxM_LH; - private: - /// 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(); }; /** Mass-conservative prolongation operator going in the opposite direction From d4d1eadb1ab7f73d28004b3a5774bf9e45b30594 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Wed, 9 Aug 2023 11:21:21 -0700 Subject: [PATCH 29/34] Avoid GetRestrictionTransposeOperator in L2ProjectionH1Space See the bug fixed in PR #3821 --- fem/transfer.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 4ec66ea633..d2b1ee6ad8 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -956,10 +956,10 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofsTranspose( void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofsTranspose( const FiniteElementSpace& fes, const Vector &X, Vector& x) const { - const Operator* Rtranspose = fes.GetRestrictionTransposeOperator(); - if (Rtranspose) + const Operator *R = fes.GetRestrictionOperator(); + if (R) { - Rtranspose->Mult(X, x); + R->MultTranspose(X, x); } else { From 472f2b83d2cd694e20ad0b6bc88ed9b938d2a836 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Wed, 9 Aug 2023 11:22:45 -0700 Subject: [PATCH 30/34] Small changes to dual vector transfer in lor-transfer miniapps --- miniapps/tools/lor-transfer.cpp | 13 +++++----- miniapps/tools/plor-transfer.cpp | 44 +++++++++++++++----------------- 2 files changed, 27 insertions(+), 30 deletions(-) diff --git a/miniapps/tools/lor-transfer.cpp b/miniapps/tools/lor-transfer.cpp index 7c2269a411..607768b445 100644 --- a/miniapps/tools/lor-transfer.cpp +++ b/miniapps/tools/lor-transfer.cpp @@ -151,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; } @@ -187,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 @@ -233,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; diff --git a/miniapps/tools/plor-transfer.cpp b/miniapps/tools/plor-transfer.cpp index b0c3d00d41..d8d87b8e86 100644 --- a/miniapps/tools/plor-transfer.cpp +++ b/miniapps/tools/plor-transfer.cpp @@ -162,6 +162,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; } @@ -219,22 +223,16 @@ int main(int argc, char *argv[]) }; 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(); - Vector rho_true(rho.ParFESpace()->GetTrueVSize()); - rho.GetTrueDofs(rho_true); - Vector M_rho_true(M_rho.ParFESpace()->GetTrueVSize()); - M_ho_tdof->Mult(rho_true, M_rho_true); - M_rho.ParFESpace()->GetRestrictionTransposeOperator()->Mult(M_rho_true, M_rho); P.MultTranspose(M_rho, M_rho_lor); - Vector M_rho_lor_true(M_rho_lor.ParFESpace()->GetTrueVSize()); - M_rho_lor.ParFESpace()->GetProlongationMatrix()->MultTranspose(M_rho_lor, - M_rho_lor_true); - double ho_dual_mass = global_sum(M_rho_true); - double lor_dual_mass = global_sum(M_rho_lor_true); + 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: " << fabs(ho_dual_mass - lor_dual_mass) << endl - << endl; + cout << "HO -> LOR dual field: " << abs(ho_dual_mass - lor_dual_mass) << "\n\n"; } } @@ -275,20 +273,20 @@ int main(int argc, char *argv[]) // LOR* to HO* dual fields if (!use_pointwise_transfer) { - Vector rho_lor_true(rho_lor.ParFESpace()->GetTrueVSize()); - rho_lor.GetTrueDofs(rho_lor_true); - Vector M_rho_lor_true(M_rho_lor.ParFESpace()->GetTrueVSize()); - M_lor_tdof->Mult(rho_lor_true, M_rho_lor_true); - M_rho_lor.ParFESpace()->GetRestrictionTransposeOperator()->Mult(M_rho_lor_true, - M_rho_lor); + 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); - Vector M_rho_true(M_rho.ParFESpace()->GetTrueVSize()); - M_rho.ParFESpace()->GetProlongationMatrix()->MultTranspose(M_rho, M_rho_true); - double ho_dual_mass = global_sum(M_rho_true); - double lor_dual_mass = global_sum(M_rho_lor_true); + 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: " << fabs(ho_dual_mass - lor_dual_mass) << '\n'; + cout << "LOR -> HO dual field: " << abs(ho_dual_mass - lor_dual_mass) << '\n'; } } From c1779a49cb0cf69c3f9b8d4e1e82c2265b6c10ea Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Wed, 9 Aug 2023 11:37:17 -0700 Subject: [PATCH 31/34] Fix shadow warning --- fem/transfer.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index d2b1ee6ad8..9ca9e5f3cc 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -956,10 +956,10 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofsTranspose( void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofsTranspose( const FiniteElementSpace& fes, const Vector &X, Vector& x) const { - const Operator *R = fes.GetRestrictionOperator(); - if (R) + const Operator *R_op = fes.GetRestrictionOperator(); + if (R_op) { - R->MultTranspose(X, x); + R_op->MultTranspose(X, x); } else { From 8806eec449e5a72985f362278f6a77264745c398 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Wed, 9 Aug 2023 11:55:31 -0700 Subject: [PATCH 32/34] add doxygen comments --- fem/transfer.cpp | 12 ++++++------ fem/transfer.hpp | 27 ++++++++++++++++++++------- 2 files changed, 26 insertions(+), 13 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 9ca9e5f3cc..1a2fc8156e 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -940,30 +940,30 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofs( } void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofsTranspose( - const FiniteElementSpace& fes, const Vector& x, Vector& X) const + const FiniteElementSpace& fes, const Vector& b, Vector& B) const { const Operator* P = fes.GetProlongationMatrix(); if (P) { - P->MultTranspose(x, X); + P->MultTranspose(b, B); } else { - X = x; + B = b; } } void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofsTranspose( - const FiniteElementSpace& fes, const Vector &X, Vector& x) const + const FiniteElementSpace& fes, const Vector &B, Vector& b) const { const Operator *R_op = fes.GetRestrictionOperator(); if (R_op) { - R_op->MultTranspose(X, x); + R_op->MultTranspose(B, b); } else { - x = X; + b = B; } } diff --git a/fem/transfer.hpp b/fem/transfer.hpp index 2745e8b092..ddbf1d5c0c 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -299,25 +299,38 @@ protected: virtual void SetRelTol(double p_rtol_); virtual void SetAbsTol(double p_atol_); protected: - /// Sets up the PCG solver (sets parameters, operator, and preconditioner) + /// @brief Sets up the PCG solver (sets parameters, operator, and + /// preconditioner) void SetupPCG(); - /// Computes on-rank R and M_LH matrices + /// @brief 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; + /// @brief 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; + const Vector& b, + Vector& B) 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; + const Vector& B, + Vector& b) 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 + /// @brief Returns the inverse of an on-rank lumped mass matrix + /// void LumpedMassInverse(Vector& ML_inv) const; /// Computes sparsity pattern and initializes R matrix. Based on /// BilinearForm::AllocMat() except maps between coarse HO elements and From f986b2022bca10f2f66931ce3cec2e06fcdc3d76 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Wed, 9 Aug 2023 12:06:29 -0700 Subject: [PATCH 33/34] Minor Doxygen formatting --- fem/transfer.hpp | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/fem/transfer.hpp b/fem/transfer.hpp index ddbf1d5c0c..caea4e78f9 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,8 +255,8 @@ 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. }; /** Projection operator between a H1 high-order finite element space on a @@ -299,18 +305,15 @@ protected: virtual void SetRelTol(double p_rtol_); virtual void SetAbsTol(double p_atol_); protected: - /// @brief Sets up the PCG solver (sets parameters, operator, and - /// preconditioner) + /// Sets up the PCG solver (sets parameters, operator, and preconditioner) void SetupPCG(); - /// @brief Computes on-rank R and M_LH matrices - /// + /// 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; - /// @brief Sets dof values given a vector of tdofs and a finite element - /// space + /// Sets dof values given a vector of tdofs and a finite element space void SetFromTDofs(const FiniteElementSpace& fes, const Vector& X, Vector& x) const; @@ -329,12 +332,12 @@ protected: void TDofsListByVDim(const FiniteElementSpace& fes, int vdim, Array& vdofs_list) const; - /// @brief Returns the inverse of an on-rank lumped mass matrix - /// + /// Returns the inverse of an on-rank lumped mass matrix void LumpedMassInverse(Vector& ML_inv) const; - /// Computes sparsity pattern and initializes R matrix. Based on - /// BilinearForm::AllocMat() except maps between coarse HO elements and - /// refined LOR elements. + /// @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; From 94fb0f94c51c7973f18638ce426ebe13631414b5 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Wed, 9 Aug 2023 12:34:22 -0700 Subject: [PATCH 34/34] Fix shadow warning --- fem/transfer.cpp | 12 ++++++------ fem/transfer.hpp | 8 ++++---- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 1a2fc8156e..9ca9e5f3cc 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -940,30 +940,30 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofs( } void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofsTranspose( - const FiniteElementSpace& fes, const Vector& b, Vector& B) const + const FiniteElementSpace& fes, const Vector& x, Vector& X) const { const Operator* P = fes.GetProlongationMatrix(); if (P) { - P->MultTranspose(b, B); + P->MultTranspose(x, X); } else { - B = b; + X = x; } } void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofsTranspose( - const FiniteElementSpace& fes, const Vector &B, Vector& b) const + const FiniteElementSpace& fes, const Vector &X, Vector& x) const { const Operator *R_op = fes.GetRestrictionOperator(); if (R_op) { - R_op->MultTranspose(B, b); + R_op->MultTranspose(X, x); } else { - b = B; + x = X; } } diff --git a/fem/transfer.hpp b/fem/transfer.hpp index caea4e78f9..553f037cc6 100644 --- a/fem/transfer.hpp +++ b/fem/transfer.hpp @@ -320,13 +320,13 @@ protected: /// @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& b, - Vector& B) const; + 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& B, - Vector& b) const; + 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,