From 8bc17f0f8109a3a6dd7e9fb60c3415782ca2535c Mon Sep 17 00:00:00 2001 From: Christian Glusa Date: Thu, 5 Dec 2024 17:39:23 -0700 Subject: [PATCH] MueLu RefMaxwell: Add option to pass in material Signed-off-by: Christian Glusa --- .../src/Operators/MueLu_RefMaxwell_decl.hpp | 60 +++++++++++++++++-- .../src/Operators/MueLu_RefMaxwell_def.hpp | 29 ++++++--- packages/muelu/test/maxwell/Maxwell3D.cpp | 2 +- 3 files changed, 77 insertions(+), 14 deletions(-) diff --git a/packages/muelu/src/Operators/MueLu_RefMaxwell_decl.hpp b/packages/muelu/src/Operators/MueLu_RefMaxwell_decl.hpp index bead022651eb..ff143ef7459d 100644 --- a/packages/muelu/src/Operators/MueLu_RefMaxwell_decl.hpp +++ b/packages/muelu/src/Operators/MueLu_RefMaxwell_decl.hpp @@ -264,6 +264,36 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator &SM_Matrix, + const Teuchos::RCP &Dk_1, + const Teuchos::RCP &Dk_2, + const Teuchos::RCP &D0, + const Teuchos::RCP &M1_beta, + const Teuchos::RCP &M1_alpha, + const Teuchos::RCP &Mk_one, + const Teuchos::RCP &Mk_1_one, + const Teuchos::RCP &invMk_1_invBeta, + const Teuchos::RCP &invMk_2_invAlpha, + const Teuchos::RCP &Nullspace11, + const Teuchos::RCP &Nullspace22, + const Teuchos::RCP &NodalCoords, + const Teuchos::RCP &Material_beta, + const Teuchos::RCP &Material_alpha, + Teuchos::ParameterList &List, + bool ComputePrec = true) { + int spaceNumber = List.get("refmaxwell: space number", 1); + initialize(spaceNumber, + Dk_1, Dk_2, D0, + M1_beta, M1_alpha, + Mk_one, Mk_1_one, + invMk_1_invBeta, invMk_2_invAlpha, + Nullspace11, Nullspace22, NodalCoords, + Material_beta, Material_alpha, List); resetMatrix(SM_Matrix, ComputePrec); } @@ -289,7 +319,21 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator &NodalCoords, Teuchos::ParameterList &List, bool ComputePrec = true) { - initialize(D0_Matrix, Ms_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, List); + initialize(D0_Matrix, Ms_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, Teuchos::null, List); + resetMatrix(SM_Matrix, ComputePrec); + } + + RefMaxwell(const Teuchos::RCP &SM_Matrix, + const Teuchos::RCP &D0_Matrix, + const Teuchos::RCP &Ms_Matrix, + const Teuchos::RCP &M0inv_Matrix, + const Teuchos::RCP &M1_Matrix, + const Teuchos::RCP &Nullspace11, + const Teuchos::RCP &NodalCoords, + const Teuchos::RCP &Material, + Teuchos::ParameterList &List, + bool ComputePrec = true) { + initialize(D0_Matrix, Ms_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, Material, List); resetMatrix(SM_Matrix, ComputePrec); } @@ -312,7 +356,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator &NodalCoords, Teuchos::ParameterList &List, bool ComputePrec = true) { - initialize(D0_Matrix, M1_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, List); + initialize(D0_Matrix, M1_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, Teuchos::null, List); resetMatrix(SM_Matrix, ComputePrec); } @@ -332,7 +376,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator &NodalCoords, Teuchos::ParameterList &List) : SM_Matrix_(Teuchos::null) { - initialize(D0_Matrix, M1_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, List); + initialize(D0_Matrix, M1_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, Teuchos::null, List); } /** Constructor with Jacobian (no add on) @@ -352,7 +396,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator &NodalCoords, Teuchos::ParameterList &List, bool ComputePrec) { - initialize(D0_Matrix, M1_Matrix, Teuchos::null, M1_Matrix, Nullspace11, NodalCoords, List); + initialize(D0_Matrix, M1_Matrix, Teuchos::null, M1_Matrix, Nullspace11, NodalCoords, Teuchos::null, List); resetMatrix(SM_Matrix, ComputePrec); } @@ -370,7 +414,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator &NodalCoords, Teuchos::ParameterList &List) : SM_Matrix_(Teuchos::null) { - initialize(D0_Matrix, M1_Matrix, Teuchos::null, M1_Matrix, Nullspace11, NodalCoords, List); + initialize(D0_Matrix, M1_Matrix, Teuchos::null, M1_Matrix, Nullspace11, NodalCoords, Teuchos::null, List); } /** Constructor with parameter list @@ -449,6 +493,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator &M1_Matrix, const Teuchos::RCP &Nullspace11, const Teuchos::RCP &NodalCoords, + const Teuchos::RCP &Material, Teuchos::ParameterList &List); /** Initialize with matrices except the Jacobian (don't compute the preconditioner) @@ -481,6 +526,8 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator &Nullspace11, const Teuchos::RCP &Nullspace22, const Teuchos::RCP &NodalCoords, + const Teuchos::RCP &Material_beta, + const Teuchos::RCP &Material_alpha, Teuchos::ParameterList &List); //! Determine how large the sub-communicators for the two hierarchies should be @@ -549,6 +596,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator &A, const Teuchos::RCP &Nullspace, const Teuchos::RCP &Coords, + const Teuchos::RCP &Material, Teuchos::ParameterList ¶ms, std::string &label, const bool reuse, @@ -612,6 +660,8 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator Mk_one_, Mk_1_one_; //! mass matrices on first space with weights beta and alpha respectively Teuchos::RCP M1_beta_, M1_alpha_; + //! material for first space + Teuchos::RCP Material_beta_, Material_alpha_; //! special prolongator for 11 block and its transpose Teuchos::RCP P11_, R11_; //! special prolongator for 22 block and its transpose diff --git a/packages/muelu/src/Operators/MueLu_RefMaxwell_def.hpp b/packages/muelu/src/Operators/MueLu_RefMaxwell_def.hpp index 54ba2b31c9d7..247d00743533 100644 --- a/packages/muelu/src/Operators/MueLu_RefMaxwell_def.hpp +++ b/packages/muelu/src/Operators/MueLu_RefMaxwell_def.hpp @@ -286,11 +286,11 @@ void RefMaxwell::setParameters(Teucho // This should be taken out again as soon as // CoalesceDropFactory_kokkos supports BlockSize > 1 and // drop tol != 0.0 - if (useKokkos_ && precList11_.isParameter("aggregation: drop tol") && precList11_.get("aggregation: drop tol") != 0.0) { - GetOStream(Warnings0) << solverName_ + "::compute(): Setting \"aggregation: drop tol\". to 0.0, since CoalesceDropFactory_kokkos does not " - << "support BlockSize > 1 and drop tol != 0.0" << std::endl; - precList11_.set("aggregation: drop tol", 0.0); - } + // if (useKokkos_ && precList11_.isParameter("aggregation: drop tol") && precList11_.get("aggregation: drop tol") != 0.0) { + // GetOStream(Warnings0) << solverName_ + "::compute(): Setting \"aggregation: drop tol\". to 0.0, since CoalesceDropFactory_kokkos does not " + // << "support BlockSize > 1 and drop tol != 0.0" << std::endl; + // precList11_.set("aggregation: drop tol", 0.0); + // } } template @@ -469,7 +469,7 @@ void RefMaxwell::compute(bool reuse) if (!coarseA11_.is_null()) { VerbLevel verbosityLevel = VerboseObject::GetDefaultVerbLevel(); std::string label("coarseA11"); - setupSubSolve(HierarchyCoarse11_, thyraPrecOpH_, coarseA11_, NullspaceCoarse11_, CoordsCoarse11_, precList11_, label, reuse); + setupSubSolve(HierarchyCoarse11_, thyraPrecOpH_, coarseA11_, NullspaceCoarse11_, CoordsCoarse11_, Material_beta_, precList11_, label, reuse); VerboseObject::SetDefaultVerbLevel(verbosityLevel); } @@ -525,9 +525,9 @@ void RefMaxwell::compute(bool reuse) int maxLevels = precList22_.get("max levels", MasterList::getDefault("max levels")); if (maxLevels < 2) precList22_.set("max levels", 2); - setupSubSolve(Hierarchy22_, thyraPrecOp22_, A22_, Teuchos::null, Teuchos::null, precList22_, label, reuse, /*isSingular=*/globalNumberBoundaryUnknowns11_ == 0); + setupSubSolve(Hierarchy22_, thyraPrecOp22_, A22_, Teuchos::null, Teuchos::null, Material_alpha_, precList22_, label, reuse, /*isSingular=*/globalNumberBoundaryUnknowns11_ == 0); } else - setupSubSolve(Hierarchy22_, thyraPrecOp22_, A22_, CoarseNullspace22_, Coords22_, precList22_, label, reuse, /*isSingular=*/globalNumberBoundaryUnknowns11_ == 0); + setupSubSolve(Hierarchy22_, thyraPrecOp22_, A22_, CoarseNullspace22_, Coords22_, Material_alpha_, precList22_, label, reuse, /*isSingular=*/globalNumberBoundaryUnknowns11_ == 0); VerboseObject::SetDefaultVerbLevel(verbosityLevel); } @@ -1994,6 +1994,7 @@ void RefMaxwell::setupSubSolve(Teucho const Teuchos::RCP &A, const Teuchos::RCP &Nullspace, const Teuchos::RCP &Coords, + const Teuchos::RCP &Material, Teuchos::ParameterList ¶ms, std::string &label, const bool reuse, @@ -2013,6 +2014,8 @@ void RefMaxwell::setupSubSolve(Teucho ParameterList &userParamList = params.sublist("Preconditioner Types").sublist("MueLu").sublist("user data"); if (!Nullspace.is_null()) userParamList.set >("Nullspace", Nullspace); + if (!Material.is_null()) + userParamList.set >("Material", Material); userParamList.set >("Coordinates", Coords); } thyraPrecOp = rcp(new XpetraThyraLinearOp(coarseA11_, rcp(¶ms, false))); @@ -2027,6 +2030,8 @@ void RefMaxwell::setupSubSolve(Teucho userParamList.set >("Coordinates", Coords); if (!Nullspace.is_null()) userParamList.set >("Nullspace", Nullspace); + if (!Material.is_null()) + userParamList.set >("Material", Material); if (isSingular) { std::string coarseType = ""; @@ -2441,6 +2446,7 @@ RefMaxwell:: invMk_1_invBeta, invMk_2_invAlpha, Nullspace11, Nullspace22, NodalCoords, + Teuchos::null, Teuchos::null, List); if (SM_Matrix != Teuchos::null) @@ -2455,6 +2461,7 @@ void RefMaxwell:: const Teuchos::RCP &M1_Matrix, const Teuchos::RCP &Nullspace11, const Teuchos::RCP &NodalCoords, + const Teuchos::RCP &Material, Teuchos::ParameterList &List) { initialize(1, D0_Matrix, Teuchos::null, D0_Matrix, @@ -2463,6 +2470,7 @@ void RefMaxwell:: M0inv_Matrix, Teuchos::null, Nullspace11, Teuchos::null, NodalCoords, + Teuchos::null, Material, List); } @@ -2481,6 +2489,8 @@ void RefMaxwell:: const Teuchos::RCP &Nullspace11, const Teuchos::RCP &Nullspace22, const Teuchos::RCP &NodalCoords, + const Teuchos::RCP &Material_beta, + const Teuchos::RCP &Material_alpha, Teuchos::ParameterList &List) { spaceNumber_ = k; if (spaceNumber_ == 1) @@ -2645,6 +2655,9 @@ void RefMaxwell:: M1_beta_ = M1_beta; M1_alpha_ = M1_alpha; + Material_beta_ = Material_beta; + Material_alpha_ = Material_alpha; + Mk_one_ = Mk_one; Mk_1_one_ = Mk_1_one; diff --git a/packages/muelu/test/maxwell/Maxwell3D.cpp b/packages/muelu/test/maxwell/Maxwell3D.cpp index 283c55baec8a..7bb3bcc207d0 100644 --- a/packages/muelu/test/maxwell/Maxwell3D.cpp +++ b/packages/muelu/test/maxwell/Maxwell3D.cpp @@ -254,7 +254,7 @@ bool SetupSolve(std::map inputs) { RCP preconditioner; if (precType == "MueLu-RefMaxwell") { preconditioner = rcp(new MueLu::RefMaxwell(SM_Matrix, D0_Matrix, Ms_Matrix, M0inv_Matrix, - M1_Matrix, nullspace, coords, params)); + 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(SM_Matrix, D0_Matrix, Kn_Matrix, nullspace, coords, params));