Skip to content

Commit

Permalink
MueLu RefMaxwell: Add option to pass in material
Browse files Browse the repository at this point in the history 
Signed-off-by: Christian Glusa <[email protected]>
  • Loading branch information
@cgcgcg
cgcgcg committed Dec 10, 2024
1 parent 77e6090 commit 8bc17f0
Show file tree
Hide file tree
Showing 3 changed files with 77 additions and 14 deletions.
60 changes: 55 additions & 5 deletions packages/muelu/src/Operators/MueLu_RefMaxwell_decl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,6 +264,36 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOr
Mk_one, Mk_1_one,
invMk_1_invBeta, invMk_2_invAlpha,
Nullspace11, Nullspace22, NodalCoords,
Teuchos::null, Teuchos::null,
List);
resetMatrix(SM_Matrix, ComputePrec);
}

RefMaxwell(const Teuchos::RCP<Matrix> &SM_Matrix,
const Teuchos::RCP<Matrix> &Dk_1,
const Teuchos::RCP<Matrix> &Dk_2,
const Teuchos::RCP<Matrix> &D0,
const Teuchos::RCP<Matrix> &M1_beta,
const Teuchos::RCP<Matrix> &M1_alpha,
const Teuchos::RCP<Matrix> &Mk_one,
const Teuchos::RCP<Matrix> &Mk_1_one,
const Teuchos::RCP<Matrix> &invMk_1_invBeta,
const Teuchos::RCP<Matrix> &invMk_2_invAlpha,
const Teuchos::RCP<MultiVector> &Nullspace11,
const Teuchos::RCP<MultiVector> &Nullspace22,
const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
const Teuchos::RCP<MultiVector> &Material_beta,
const Teuchos::RCP<MultiVector> &Material_alpha,
Teuchos::ParameterList &List,
bool ComputePrec = true) {
int spaceNumber = List.get<int>("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);
}
Expand All @@ -289,7 +319,21 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOr
const Teuchos::RCP<RealValuedMultiVector> &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<Matrix> &SM_Matrix,
const Teuchos::RCP<Matrix> &D0_Matrix,
const Teuchos::RCP<Matrix> &Ms_Matrix,
const Teuchos::RCP<Matrix> &M0inv_Matrix,
const Teuchos::RCP<Matrix> &M1_Matrix,
const Teuchos::RCP<MultiVector> &Nullspace11,
const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
const Teuchos::RCP<MultiVector> &Material,
Teuchos::ParameterList &List,
bool ComputePrec = true) {
initialize(D0_Matrix, Ms_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, Material, List);
resetMatrix(SM_Matrix, ComputePrec);
}

Expand All @@ -312,7 +356,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOr
const Teuchos::RCP<RealValuedMultiVector> &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);
}

Expand All @@ -332,7 +376,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOr
const Teuchos::RCP<RealValuedMultiVector> &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)
Expand All @@ -352,7 +396,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOr
const Teuchos::RCP<RealValuedMultiVector> &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);
}

Expand All @@ -370,7 +414,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOr
const Teuchos::RCP<RealValuedMultiVector> &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
Expand Down Expand Up @@ -449,6 +493,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOr
const Teuchos::RCP<Matrix> &M1_Matrix,
const Teuchos::RCP<MultiVector> &Nullspace11,
const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
const Teuchos::RCP<MultiVector> &Material,
Teuchos::ParameterList &List);

/** Initialize with matrices except the Jacobian (don't compute the preconditioner)
Expand Down Expand Up @@ -481,6 +526,8 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOr
const Teuchos::RCP<MultiVector> &Nullspace11,
const Teuchos::RCP<MultiVector> &Nullspace22,
const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
const Teuchos::RCP<MultiVector> &Material_beta,
const Teuchos::RCP<MultiVector> &Material_alpha,
Teuchos::ParameterList &List);

//! Determine how large the sub-communicators for the two hierarchies should be
Expand Down Expand Up @@ -549,6 +596,7 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOr
const Teuchos::RCP<Matrix> &A,
const Teuchos::RCP<MultiVector> &Nullspace,
const Teuchos::RCP<RealValuedMultiVector> &Coords,
const Teuchos::RCP<MultiVector> &Material,
Teuchos::ParameterList &params,
std::string &label,
const bool reuse,
Expand Down Expand Up @@ -612,6 +660,8 @@ class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOr
Teuchos::RCP<Matrix> Mk_one_, Mk_1_one_;
//! mass matrices on first space with weights beta and alpha respectively
Teuchos::RCP<Matrix> M1_beta_, M1_alpha_;
//! material for first space
Teuchos::RCP<MultiVector> Material_beta_, Material_alpha_;
//! special prolongator for 11 block and its transpose
Teuchos::RCP<Matrix> P11_, R11_;
//! special prolongator for 22 block and its transpose
Expand Down
29 changes: 21 additions & 8 deletions packages/muelu/src/Operators/MueLu_RefMaxwell_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -286,11 +286,11 @@ void RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::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<double>("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<double>("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 <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Expand Down Expand Up @@ -469,7 +469,7 @@ void RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::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);
}

Expand Down Expand Up @@ -525,9 +525,9 @@ void RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::compute(bool reuse)
int maxLevels = precList22_.get("max levels", MasterList::getDefault<int>("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);
}
Expand Down Expand Up @@ -1994,6 +1994,7 @@ void RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupSubSolve(Teucho
const Teuchos::RCP<Matrix> &A,
const Teuchos::RCP<MultiVector> &Nullspace,
const Teuchos::RCP<RealValuedMultiVector> &Coords,
const Teuchos::RCP<MultiVector> &Material,
Teuchos::ParameterList &params,
std::string &label,
const bool reuse,
Expand All @@ -2013,6 +2014,8 @@ void RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupSubSolve(Teucho
ParameterList &userParamList = params.sublist("Preconditioner Types").sublist("MueLu").sublist("user data");
if (!Nullspace.is_null())
userParamList.set<RCP<MultiVector> >("Nullspace", Nullspace);
if (!Material.is_null())
userParamList.set<RCP<MultiVector> >("Material", Material);
userParamList.set<RCP<RealValuedMultiVector> >("Coordinates", Coords);
}
thyraPrecOp = rcp(new XpetraThyraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(coarseA11_, rcp(&params, false)));
Expand All @@ -2027,6 +2030,8 @@ void RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupSubSolve(Teucho
userParamList.set<RCP<RealValuedMultiVector> >("Coordinates", Coords);
if (!Nullspace.is_null())
userParamList.set<RCP<MultiVector> >("Nullspace", Nullspace);
if (!Material.is_null())
userParamList.set<RCP<MultiVector> >("Material", Material);

if (isSingular) {
std::string coarseType = "";
Expand Down Expand Up @@ -2441,6 +2446,7 @@ RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
invMk_1_invBeta, invMk_2_invAlpha,
Nullspace11, Nullspace22,
NodalCoords,
Teuchos::null, Teuchos::null,
List);

if (SM_Matrix != Teuchos::null)
Expand All @@ -2455,6 +2461,7 @@ void RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
const Teuchos::RCP<Matrix> &M1_Matrix,
const Teuchos::RCP<MultiVector> &Nullspace11,
const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
const Teuchos::RCP<MultiVector> &Material,
Teuchos::ParameterList &List) {
initialize(1,
D0_Matrix, Teuchos::null, D0_Matrix,
Expand All @@ -2463,6 +2470,7 @@ void RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
M0inv_Matrix, Teuchos::null,
Nullspace11, Teuchos::null,
NodalCoords,
Teuchos::null, Material,
List);
}

Expand All @@ -2481,6 +2489,8 @@ void RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
const Teuchos::RCP<MultiVector> &Nullspace11,
const Teuchos::RCP<MultiVector> &Nullspace22,
const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
const Teuchos::RCP<MultiVector> &Material_beta,
const Teuchos::RCP<MultiVector> &Material_alpha,
Teuchos::ParameterList &List) {
spaceNumber_ = k;
if (spaceNumber_ == 1)
Expand Down Expand Up @@ -2645,6 +2655,9 @@ void RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
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;

Expand Down
2 changes: 1 addition & 1 deletion packages/muelu/test/maxwell/Maxwell3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,7 +254,7 @@ bool SetupSolve(std::map<std::string, void*> inputs) {
RCP<Operator> preconditioner;
if (precType == "MueLu-RefMaxwell") {
preconditioner = rcp(new MueLu::RefMaxwell<SC, LO, GO, NO>(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<SC, LO, GO, NO>(SM_Matrix, D0_Matrix, Kn_Matrix, nullspace, coords, params));
Expand Down

0 comments on commit 8bc17f0

Please sign in to comment.