Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #ifdef MOOSE_MFEM_ENABLED
11 :
12 : #include "MFEMGeometricMultigridSolver.h"
13 : #include "MFEMProblem.h"
14 : #include "EquationSystem.h"
15 :
16 : registerMooseObject("MooseApp", MFEMGeometricMultigridSolver);
17 :
18 9 : MFEMGeometricMultigridSolver::MGProxy::MGProxy(MFEMGeometricMultigridSolver & owner) : _owner(owner)
19 : {
20 : // MGProxy is installed as a preconditioner, so Mult() should overwrite its output vector rather
21 : // than treating it as an initial iterate. This matches MFEM's
22 : // IterativeSolver::SetPreconditioner() convention, which sets the preconditioner's iterative_mode
23 : // to false.
24 9 : iterative_mode = false;
25 9 : }
26 :
27 : void
28 7 : MFEMGeometricMultigridSolver::MGProxy::SetMG(mfem::GeometricMultigrid & mg)
29 : {
30 7 : _mg = &mg;
31 7 : height = mg.Height();
32 7 : width = mg.Width();
33 7 : }
34 :
35 : void
36 9 : MFEMGeometricMultigridSolver::MGProxy::SetOperator(const mfem::Operator & op)
37 : {
38 9 : _owner.BuildMultigrid(op);
39 7 : }
40 :
41 : void
42 35 : MFEMGeometricMultigridSolver::MGProxy::Mult(const mfem::Vector & x, mfem::Vector & y) const
43 : {
44 35 : MFEM_VERIFY(_mg, "MGProxy: GeometricMultigrid not yet built");
45 35 : _mg->Mult(x, y);
46 35 : }
47 :
48 : InputParameters
49 2146 : MFEMGeometricMultigridSolver::validParams()
50 : {
51 2146 : InputParameters params = Moose::MFEM::LinearSolverBase::validParams();
52 4292 : params.addClassDescription(
53 : "Geometric (p-)multigrid preconditioner backed by mfem::GeometricMultigrid. "
54 : "Requires a linear equation system, an MFEMFESpaceHierarchy, and per-level smoother "
55 : "objects.");
56 :
57 8584 : params.addRequiredParam<std::string>("variable",
58 : "Name of the trial variable this preconditioner acts on.");
59 8584 : params.addRequiredParam<std::vector<MFEMSolverName>>(
60 : "smoothers",
61 : "Names of LinearSolverBase objects used as smoothers on the interior levels "
62 : "(levels 1 to N-1). May have length 1 (used on all interior levels) or "
63 : "N-1 (one per interior level, ordered coarse-to-fine).");
64 8584 : params.addRequiredParam<MFEMSolverName>(
65 : "coarse_solver", "Name of the LinearSolverBase used on the coarsest level.");
66 8584 : params.addParam<std::vector<std::string>>(
67 : "assembly_levels",
68 : {"legacy"},
69 : "Assembly level for each level in the hierarchy. Valid values: 'legacy', 'full', "
70 : "'element', 'partial', 'none'. May have length 1 (used on all N levels) or N.");
71 2146 : return params;
72 6438 : }
73 :
74 11 : MFEMGeometricMultigridSolver::MFEMGeometricMultigridSolver(const InputParameters & parameters)
75 : : Moose::MFEM::LinearSolverBase(parameters),
76 11 : _var_name(getParam<std::string>("variable")),
77 22 : _smoother_names(getParam<std::vector<MFEMSolverName>>("smoothers")),
78 55 : _coarse_solver_name(getParam<MFEMSolverName>("coarse_solver"))
79 : {
80 11 : auto & problem = getMFEMProblem();
81 11 : auto eq_sys = problem.getProblemData().eqn_system;
82 :
83 11 : if (eq_sys->IsEigen() || eq_sys->IsComplex())
84 2 : mooseError("GeometricMultigridSolver '", name(), "': requires a real, non-eigen eq. system");
85 :
86 : // Co-own the hierarchy so it outlives this solver.
87 27 : if (auto * hierarchy_name = problem.getMFEMObject<MFEMVariable>("MooseVariableBase", _var_name)
88 36 : .queryParam<std::string>("fespace_hierarchy"))
89 9 : _hierarchy = problem.getProblemData().fespace_hierarchies.GetShared(*hierarchy_name);
90 : else
91 0 : paramError("variable", "must be associated with an MFEMFESpaceHierarchy.");
92 :
93 : // Parse assembly levels, optionally expanding a single input value to all levels.
94 9 : const int N = _hierarchy->GetNumLevels();
95 : mooseAssert(N, "Malformed MFEMFESpaceHierarchy w/ no levels");
96 18 : const auto & asm_strs = getParam<std::vector<std::string>>("assembly_levels");
97 9 : const int n_asm = asm_strs.size();
98 9 : if (n_asm != 1 && n_asm != N)
99 0 : paramError(
100 : "assembly_levels", "must have length 1 or N = ", N, " (total levels), got ", n_asm, ".");
101 :
102 9 : _assembly_levels.resize(N);
103 27 : for (const auto i : make_range(N))
104 18 : _assembly_levels[i] = ParseAssemblyLevel(n_asm == 1 ? asm_strs[0] : asm_strs[i]);
105 :
106 9 : ConstructSolver();
107 9 : }
108 :
109 : void
110 9 : MFEMGeometricMultigridSolver::ConstructSolver()
111 : {
112 9 : _mg.reset();
113 9 : _level_ops.clear();
114 9 : _level_blfs.clear();
115 :
116 9 : auto proxy = std::make_unique<MGProxy>(*this);
117 9 : _mg_proxy = proxy.get();
118 9 : _solver = std::move(proxy);
119 9 : }
120 :
121 : mfem::AssemblyLevel
122 18 : MFEMGeometricMultigridSolver::ParseAssemblyLevel(const std::string & s) const
123 : {
124 54 : static MooseEnum assembly_levels("legacy full element partial none", "legacy");
125 18 : return (assembly_levels = s).getEnum<mfem::AssemblyLevel>();
126 : }
127 :
128 : void
129 0 : MFEMGeometricMultigridSolver::SetOperator(mfem::Operator & op)
130 : {
131 0 : BuildMultigrid(op);
132 0 : }
133 :
134 : void
135 9 : MFEMGeometricMultigridSolver::BuildMultigrid(const mfem::Operator & op)
136 : {
137 9 : auto & problem = getMFEMProblem();
138 9 : auto eq_sys = problem.getProblemData().eqn_system;
139 :
140 9 : if (eq_sys->IsNonlinear() || eq_sys->IsMultivariate())
141 2 : mooseError("GeometricMultigridSolver '", name(), "': requires a univariate, linear eq. system");
142 :
143 7 : const int N = _hierarchy->GetNumLevels();
144 7 : const int finest_level = _hierarchy->GetFinestLevelIndex();
145 :
146 : // Validate smoother vector length (levels 1 to N-1 each need a smoother).
147 7 : const int n_smooth = _smoother_names.size();
148 7 : if (n_smooth != 1 && n_smooth != N - 1)
149 0 : paramError("smoothers", "must have length 1 or N-1 = ", N - 1, ", got ", n_smooth, ".");
150 :
151 14 : auto get_smoother = [&](int level) -> Moose::MFEM::LinearSolverBase &
152 : {
153 14 : if (level == 0)
154 28 : return problem.getMFEMObject<Moose::MFEM::LinearSolverBase>("Moose::MFEM::SolverBase",
155 7 : _coarse_solver_name);
156 7 : const std::string & sname = (n_smooth == 1) ? _smoother_names[0] : _smoother_names[level - 1];
157 21 : return problem.getMFEMObject<Moose::MFEM::LinearSolverBase>("Moose::MFEM::SolverBase", sname);
158 7 : };
159 :
160 : // Obtain essential boundary attribute markers from the equation system.
161 7 : mfem::Array<int> & ess_bdr = eq_sys->GetEssentialBoundaryMarkers(_var_name);
162 :
163 : auto & finest_fespace =
164 7 : static_cast<mfem::ParFiniteElementSpace &>(_hierarchy->GetFESpaceAtLevel(finest_level));
165 7 : const int finest_size = finest_fespace.GetTrueVSize();
166 7 : if (op.Height() != finest_size || op.Width() != finest_size)
167 0 : mooseError("GeometricMultigridSolver '",
168 0 : name(),
169 : "': incoming fine operator has size ",
170 0 : op.Height(),
171 : " x ",
172 0 : op.Width(),
173 : ", but the finest hierarchy space has true size ",
174 : finest_size,
175 : ".");
176 :
177 : // Build new levels' forms; accumulate before touching _mg / _level_*.
178 7 : std::vector<std::shared_ptr<mfem::ParBilinearForm>> new_blfs;
179 7 : std::vector<std::unique_ptr<mfem::OperatorHandle>> new_level_ops;
180 7 : new_level_ops.reserve(N - 1);
181 :
182 7 : auto mg = std::make_unique<mfem::GeometricMultigrid>(*_hierarchy, ess_bdr);
183 7 : auto * mg_ptr = mg.get();
184 :
185 21 : for (const auto level : make_range(N))
186 : {
187 : auto & level_fespace =
188 14 : static_cast<mfem::ParFiniteElementSpace &>(_hierarchy->GetFESpaceAtLevel(level));
189 :
190 : // Compute essential true DoFs for this level.
191 14 : mfem::Array<int> level_tdofs;
192 14 : level_fespace.GetEssentialTrueDofs(ess_bdr, level_tdofs);
193 :
194 : // Build level operator.
195 14 : mfem::Operator * level_op = nullptr;
196 :
197 14 : if (level == finest_level)
198 7 : level_op = const_cast<mfem::Operator *>(&op);
199 : else
200 : {
201 : auto blf =
202 7 : eq_sys->BuildBilinearFormForFESpace(_var_name, level_fespace, _assembly_levels[level]);
203 :
204 7 : auto level_op_handle = std::make_unique<mfem::OperatorHandle>();
205 7 : blf->FormSystemMatrix(level_tdofs, *level_op_handle);
206 7 : level_op = level_op_handle->Ptr();
207 7 : new_level_ops.push_back(std::move(level_op_handle));
208 7 : new_blfs.push_back(std::move(blf));
209 7 : }
210 :
211 : // Configure the smoother / coarse solver with this level's operator.
212 : // Each smoother's SetOperator() owns full initialization.
213 14 : auto & level_smoother = get_smoother(level);
214 14 : level_smoother.SetOperator(*level_op);
215 :
216 28 : mg_ptr->AddLevel(
217 14 : level_op, &level_smoother.GetSolver(), /*ownOperator=*/false, /*ownSmoother=*/false);
218 14 : }
219 :
220 : // Atomically replace:
221 : // 1. Old MG freed, dropping raw pointers into level operators.
222 : // 2. Old operator handles freed before old forms they may wrap.
223 : // 3. Proxy updated to point at the new MG and level data.
224 7 : _mg = std::move(mg);
225 7 : _level_ops = std::move(new_level_ops);
226 7 : _level_blfs = std::move(new_blfs);
227 7 : _mg_proxy->SetMG(*_mg);
228 7 : }
229 : #endif
|
