Merge pull request #4271 from mfem/symmatcoeff-project-fix

Fix `SymmetricMatrixCoefficient::ProjectSymmetric` bug

.gitlab/lassen-build-and-test.yml

@@ -14,14 +14,14 @@ stages:
   - build_and_test
   - report
 
-opt_mpi_cuda_xl_16_1_1_12:
+opt_mpi_cuda_gcc:
   variables:
-    SPEC: "%xl@16.1.1.12 +mpi +cuda cuda_arch=70"
+    SPEC: "%gcc@8.3.1 +mpi +cuda cuda_arch=70"
   extends: .build_and_test_on_lassen
 
-opt_mpi_cuda_hypre_cuda_xl:
+opt_mpi_cuda_hypre_cuda_gcc:
   variables:
-    SPEC: "%xl@16.1.1.12 +mpi +cuda cuda_arch=70 ^hypre+cuda~shared cuda_arch=70"
+    SPEC: "%gcc@8.3.1 +mpi +cuda cuda_arch=70 ^hypre+cuda~shared cuda_arch=70"
   extends: .build_and_test_on_lassen
 
 # Jobs report

fem/coefficient.cpp

@@ -807,6 +807,7 @@ void SymmetricMatrixCoefficient::ProjectSymmetric(QuadratureFunction &qf)
 
    QuadratureSpaceBase &qspace = *qf.GetSpace();
    const int ne = qspace.GetNE();
+   qf.HostWrite();
    DenseMatrix values;
    DenseSymmetricMatrix matrix;
    for (int iel = 0; iel < ne; ++iel)
@@ -818,7 +819,7 @@ void SymmetricMatrixCoefficient::ProjectSymmetric(QuadratureFunction &qf)
       {
          const IntegrationPoint &ip = ir[iq];
          T.SetIntPoint(&ip);
-         matrix.UseExternalData(&values(0, iq), vdim);
+         matrix.UseExternalData(&values(0, iq), height);
          Eval(matrix, T, ip);
       }
    }
@@ -828,13 +829,12 @@ void SymmetricMatrixCoefficient::ProjectSymmetric(QuadratureFunction &qf)
 void SymmetricMatrixCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
                                       const IntegrationPoint &ip)
 {
-   mat.SetSize(height);
-   Eval(mat, T, ip);
+   Eval(mat_aux, T, ip);
    for (int j = 0; j < width; ++j)
    {
       for (int i = 0; i < height; ++ i)
       {
-         K(i, j) = mat(i, j);
+         K(i, j) = mat_aux(i, j);
       }
    }
 }

fem/coefficient.hpp

@@ -1426,12 +1426,13 @@ class SumCoefficient : public Coefficient
 class SymmetricMatrixCoefficient : public MatrixCoefficient
 {
 protected:
+
    /// Internal matrix used when evaluating this coefficient as a DenseMatrix.
-   DenseSymmetricMatrix mat;
+   mutable DenseSymmetricMatrix mat_aux;
 public:
    /// Construct a dim x dim matrix coefficient.
    explicit SymmetricMatrixCoefficient(int dimension)
-      : MatrixCoefficient(dimension, true) { }
+      : MatrixCoefficient(dimension, true), mat_aux(height) { }
 
    /// Get the size of the matrix.
    int GetSize() const { return height; }
@@ -1464,8 +1465,9 @@ class SymmetricMatrixCoefficient : public MatrixCoefficient
    virtual void Eval(DenseMatrix &K, ElementTransformation &T,
                      const IntegrationPoint &ip);
 
-   /// Return a reference to the constant matrix.
-   const DenseSymmetricMatrix& GetMatrix() { return mat; }
+
+   /// @deprecated Return a reference to the internal matrix used when evaluating this coefficient as a DenseMatrix.
+   MFEM_DEPRECATED const DenseSymmetricMatrix& GetMatrix() { return mat_aux; }
 
    virtual ~SymmetricMatrixCoefficient() { }
 };
@@ -1485,6 +1487,10 @@ class SymmetricMatrixConstantCoefficient : public SymmetricMatrixCoefficient
    /// Evaluate the matrix coefficient at @a ip.
    virtual void Eval(DenseSymmetricMatrix &M, ElementTransformation &T,
                      const IntegrationPoint &ip) { M = mat; }
+
+   /// Return a reference to the constant matrix.
+   const DenseSymmetricMatrix& GetMatrix() { return mat; }
+
 };
 
 

linalg/hypre.cpp

@@ -1016,7 +1016,8 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm,
    {
       HypreReadWrite();
       hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
-      SyncBackCSR(diag, mem_diag); // update diag, if needed
+      // update diag, if needed
+      if (!own_diag_offd) { SyncBackCSR(diag, mem_diag); }
    }
 
    hypre_MatvecCommPkgCreate(A);

linalg/symmat.cpp

@@ -58,6 +58,20 @@ DenseSymmetricMatrix &DenseSymmetricMatrix::operator=(real_t c)
    return *this;
 }
 
+DenseSymmetricMatrix &DenseSymmetricMatrix::operator=(const DenseSymmetricMatrix
+                                                      &m)
+{
+   SetSize(m.height);
+
+   const int hw = m.GetStoredSize();
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] = m.data[i];
+   }
+
+   return *this;
+}
+
 real_t &DenseSymmetricMatrix::Elem(int i, int j)
 {
    return (*this)(i,j);
@@ -89,11 +103,6 @@ MatrixInverse *DenseSymmetricMatrix::Inverse() const
    return nullptr;
 }
 
-void DenseSymmetricMatrix::Print (std::ostream & os, int width_) const
-{
-   mfem_error("DenseSymmetricMatrix::Print() not implemented!");
-}
-
 DenseSymmetricMatrix::~DenseSymmetricMatrix()
 {
    data.Delete();

linalg/symmat.hpp

@@ -102,6 +102,9 @@ class DenseSymmetricMatrix : public Matrix
 
    DenseSymmetricMatrix &operator*=(real_t c);
 
+   /// Sets the matrix size and elements equal to those of m
+   DenseSymmetricMatrix &operator=(const DenseSymmetricMatrix &m);
+
    std::size_t MemoryUsage() const { return data.Capacity() * sizeof(real_t); }
 
    /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
@@ -134,9 +137,6 @@ class DenseSymmetricMatrix : public Matrix
    /// Returns a pointer to (an approximation) of the matrix inverse.
    virtual MatrixInverse *Inverse() const;
 
-   /// Prints matrix to stream out.
-   virtual void Print (std::ostream & out = mfem::out, int width_ = 4) const;
-
    /// Destroys the symmetric matrix.
    virtual ~DenseSymmetricMatrix();
 };

miniapps/multidomain/makefile

@@ -55,6 +55,9 @@ MESH_FILES = $(notdir $(wildcard $(SRC)*.mesh))
 $(MESH_FILES): %: $(SRC)%
 	ln -sf $(<) .
 multidomain: | $(MESH_FILES)
+# The target 'copy-data' is used by the makefile in ../../tests/unit
+.PHONY: copy-data
+copy-data: | $(MESH_FILES)
 endif
 
 MFEM_TESTS = MINIAPPS

tests/gitlab/get_mfem_uberenv

@@ -18,7 +18,7 @@ set -o errexit
 set -o nounset
 
 uberenv_url="https://github.com/mfem/mfem-uberenv.git"
-uberenv_ref="67fab1adaf2d095ffa70dca97381abedbea15c89"
+uberenv_ref="fe5fa88876b29ff03177d44a5bd3e09c84ccdcbf"
 
 [[ ! -d tests/uberenv ]] && git clone ${uberenv_url} tests/uberenv
 cd tests/uberenv

tests/unit/fem/test_coefficient.cpp

@@ -303,3 +303,31 @@ TEST_CASE("Piecewise Matrix Coefficient", "[Coefficient]")
       REQUIRE(m.FNorm() == MFEM_Approx(twoNorm));
    }
 }
+
+TEST_CASE("Symmetric Matrix Coefficient", "[Coefficient]")
+{
+   int d = 3;
+   int qfdim = d*(d+1)/2;
+
+   Vector values(qfdim);
+   values.Randomize();
+
+   // Create symmetric matrix initialized w/ values
+   DenseSymmetricMatrix symMat(values.GetData(), d);
+
+   SymmetricMatrixConstantCoefficient symCoeff(symMat);
+
+   // Make mesh of size 1
+   Mesh m = Mesh::MakeCartesian1D(1);
+
+   // Define qspace on mesh w/ 1 integration point
+   QuadratureSpace qspace(&m, 1);
+
+   // Define qf
+   QuadratureFunction qf(qspace, qfdim);
+
+   symCoeff.ProjectSymmetric(qf);
+
+   // Require equality
+   REQUIRE(qf.DistanceTo(values) == MFEM_Approx(0.0));
+}

tests/unit/makefile

@@ -179,6 +179,15 @@ $(eval $(call psedov_tests,debug,DEBUG,debug))
 $(eval $(call psedov_tests,cuda,CUDA,cuda))
 $(eval $(call psedov_tests,cuda_uvm,CUDA_UVM,cuda:uvm))
 
+# For out-of-source builds, copy the meshes in ../../miniapps/multidomain from
+# the source location; these are used by 'punit_tests'.
+ifneq ($(SRC),)
+.PHONY: copy-miniapps-multidomain-data
+copy-miniapps-multidomain-data:
+	$(MAKE) -C ../../miniapps/multidomain copy-data
+punit_tests: | copy-miniapps-multidomain-data
+endif
+
 # For out-of-source builds, copy the meshes in ../../miniapps/meshing from the
 # source location; these are used by the TMOP tests.
 .PHONY: copy-miniapps-meshing-data