MueLu Maxwell1: Add option to pass in material

Signed-off-by: Christian Glusa <[email protected]>
trilinos · Dec 10, 2024 · 499d67c · 499d67c
1 parent 812f50f
commit 499d67c
Show file tree

Hide file tree

Showing 3 changed files with 35 additions and 5 deletions.
diff --git a/packages/muelu/src/Operators/MueLu_Maxwell1_decl.hpp b/packages/muelu/src/Operators/MueLu_Maxwell1_decl.hpp
@@ -84,7 +84,7 @@ class Maxwell1 : public VerboseObject, public Xpetra::Operator<Scalar, LocalOrdi
            bool ComputePrec = true)
     : mode_(MODE_STANDARD) {
     RCP<Matrix> Kn_Matrix;
-    initialize(D0_Matrix, Kn_Matrix, Nullspace, Coords, List);
+    initialize(D0_Matrix, Kn_Matrix, Nullspace, Coords, Teuchos::null, List);
     resetMatrix(SM_Matrix, ComputePrec);
   }
 
@@ -105,7 +105,29 @@ class Maxwell1 : public VerboseObject, public Xpetra::Operator<Scalar, LocalOrdi
            Teuchos::ParameterList& List,
            bool ComputePrec = true)
     : mode_(MODE_STANDARD) {
-    initialize(D0_Matrix, Kn_Matrix, Nullspace, Coords, List);
+    initialize(D0_Matrix, Kn_Matrix, Nullspace, Coords, Teuchos::null, List);
+    resetMatrix(SM_Matrix, ComputePrec);
+  }
+
+  /** Constructor with Jacobian and nodal matrix
+   *
+   * \param[in] SM_Matrix Jacobian
+   * \param[in] D0_Matrix Discrete Gradient
+   * \param[in] Kn_Matrix Nodal Laplacian
+   * \param[in] Coords Nodal coordinates
+   * \param[in] List Parameter list
+   * \param[in] ComputePrec If true, compute the preconditioner immediately
+   */
+  Maxwell1(const Teuchos::RCP<Matrix>& SM_Matrix,
+           const Teuchos::RCP<Matrix>& D0_Matrix,
+           const Teuchos::RCP<Matrix>& Kn_Matrix,
+           const Teuchos::RCP<MultiVector>& Nullspace,
+           const Teuchos::RCP<RealValuedMultiVector>& Coords,
+           const Teuchos::RCP<MultiVector>& Material,
+           Teuchos::ParameterList& List,
+           bool ComputePrec = true)
+    : mode_(MODE_STANDARD) {
+    initialize(D0_Matrix, Kn_Matrix, Nullspace, Coords, Material, List);
     resetMatrix(SM_Matrix, ComputePrec);
   }
 
@@ -127,7 +149,7 @@ class Maxwell1 : public VerboseObject, public Xpetra::Operator<Scalar, LocalOrdi
            Teuchos::ParameterList& List, const Teuchos::RCP<Matrix>& GmhdA_Matrix,
            bool ComputePrec = true)
     : mode_(MODE_GMHD_STANDARD) {
-    initialize(D0_Matrix, Kn_Matrix, Nullspace, Coords, List);
+    initialize(D0_Matrix, Kn_Matrix, Nullspace, Coords, Teuchos::null, List);
     resetMatrix(SM_Matrix, ComputePrec);
     GmhdA_Matrix_  = GmhdA_Matrix;
     HierarchyGmhd_ = rcp(new Hierarchy("HierarchyGmhd"));
@@ -151,7 +173,7 @@ class Maxwell1 : public VerboseObject, public Xpetra::Operator<Scalar, LocalOrdi
     if (List.isType<RCP<Matrix> >("Kn"))
       Kn_Matrix = List.get<RCP<Matrix> >("Kn");
 
-    initialize(D0_Matrix, Kn_Matrix, Nullspace, Coords, List);
+    initialize(D0_Matrix, Kn_Matrix, Nullspace, Coords, Teuchos::null, List);
 
     if (SM_Matrix != Teuchos::null)
       resetMatrix(SM_Matrix, ComputePrec);
@@ -215,12 +237,14 @@ class Maxwell1 : public VerboseObject, public Xpetra::Operator<Scalar, LocalOrdi
    * \param[in] Kn_Matrix Kn nodal matrix
    * \param[in] Nullspace Null space (needed for periodic)
    * \param[in] Coords Nodal coordinates
+   * \param[in] Material material multivector
    * \param[in] List Parameter list
    */
   void initialize(const Teuchos::RCP<Matrix>& D0_Matrix,
                   const Teuchos::RCP<Matrix>& Kn_Matrix,
                   const Teuchos::RCP<MultiVector>& Nullspace,
                   const Teuchos::RCP<RealValuedMultiVector>& Coords,
+                  const Teuchos::RCP<MultiVector>& Material,
                   Teuchos::ParameterList& List);
 
   //! apply RefMaxwell additive 2x2 style cycle
@@ -267,6 +291,8 @@ class Maxwell1 : public VerboseObject, public Xpetra::Operator<Scalar, LocalOrdi
   Teuchos::RCP<MultiVector> Nullspace_;
   //! Coordinates
   Teuchos::RCP<RealValuedMultiVector> Coords_;
+  //! Material
+  Teuchos::RCP<MultiVector> Material_;
   //! Some options
   bool useKokkos_, allEdgesBoundary_, allNodesBoundary_, dump_matrices_, enable_reuse_, syncTimers_;
   bool applyBCsTo22_;

diff --git a/packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp b/packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp
@@ -274,6 +274,8 @@ void Maxwell1<Scalar, LocalOrdinal, GlobalOrdinal, Node>::compute(bool reuse) {
   /* Critical ParameterList changes */
   if (!Coords_.is_null())
     precList22_.sublist("user data").set("Coordinates", Coords_);
+  if (!Material_.is_null())
+    precList22_.sublist("user data").set("Material", Material_);
 
   /* Repartitioning *must* be in sync between hierarchies, which means
      that we need to keep the importers in sync too */
@@ -766,6 +768,7 @@ void Maxwell1<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
                const Teuchos::RCP<Matrix>& Kn_Matrix,
                const Teuchos::RCP<MultiVector>& Nullspace,
                const Teuchos::RCP<RealValuedMultiVector>& Coords,
+               const Teuchos::RCP<MultiVector>& Material,
                Teuchos::ParameterList& List) {
   // some pre-conditions
   TEUCHOS_ASSERT(D0_Matrix != Teuchos::null);
@@ -827,6 +830,7 @@ void Maxwell1<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
 
   Kn_Matrix_ = Kn_Matrix;
   Coords_    = Coords;
+  Material_  = Material;
   Nullspace_ = Nullspace;
 
   if (!Kn_Matrix_.is_null()) dump(*Kn_Matrix_, "Kn.m");

diff --git a/packages/muelu/test/maxwell/Maxwell3D.cpp b/packages/muelu/test/maxwell/Maxwell3D.cpp
@@ -257,7 +257,7 @@ bool SetupSolve(std::map<std::string, void*> inputs) {
                                                                  M1_Matrix, nullspace, coords, material, params));
     } else if (precType == "MueLu-Maxwell1" || precType == "MueLu-Reitzinger") {
       if (GmhdA_Matrix.is_null())  // are we doing MHD as opposed to GMHD?
-        preconditioner = rcp(new MueLu::Maxwell1<SC, LO, GO, NO>(SM_Matrix, D0_Matrix, Kn_Matrix, nullspace, coords, params));
+        preconditioner = rcp(new MueLu::Maxwell1<SC, LO, GO, NO>(SM_Matrix, D0_Matrix, Kn_Matrix, nullspace, coords, material, params));
       else
         preconditioner = rcp(new MueLu::Maxwell1<SC, LO, GO, NO>(SM_Matrix, D0_Matrix, Kn_Matrix, nullspace, coords, params, GmhdA_Matrix));