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 2144 : MFEMGeometricMultigridSolver::validParams()
50 : {
51 2144 : InputParameters params = Moose::MFEM::LinearSolverBase::validParams();
52 4288 : 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 8576 : params.addRequiredParam<std::string>("variable",
58 : "Name of the trial variable this preconditioner acts on.");
59 8576 : params.addRequiredParam<std::string>(
60 : "fespace_hierarchy", "Name of the MFEMFESpaceHierarchy that defines the level structure.");
61 8576 : params.addRequiredParam<std::vector<MFEMSolverName>>(
62 : "smoothers",
63 : "Names of LinearSolverBase objects used as smoothers on the interior levels "
64 : "(levels 1 to N-1). May have length 1 (used on all interior levels) or "
65 : "N-1 (one per interior level, ordered coarse-to-fine).");
66 8576 : params.addRequiredParam<MFEMSolverName>(
67 : "coarse_solver", "Name of the LinearSolverBase used on the coarsest level.");
68 6432 : params.addRequiredParam<std::vector<std::string>>(
69 : "assembly_levels",
70 : "Assembly level for each level in the hierarchy. Valid values: 'legacy', 'full', "
71 : "'element', 'partial', 'none'. May have length 1 (used on all N levels) or N.");
72 2144 : return params;
73 0 : }
74 :
75 9 : MFEMGeometricMultigridSolver::MFEMGeometricMultigridSolver(const InputParameters & parameters)
76 : : Moose::MFEM::LinearSolverBase(parameters),
77 9 : _var_name(getParam<std::string>("variable")),
78 18 : _smoother_names(getParam<std::vector<MFEMSolverName>>("smoothers")),
79 45 : _coarse_solver_name(getParam<MFEMSolverName>("coarse_solver"))
80 : {
81 : // Co-own the hierarchy so it outlives this solver.
82 18 : const auto & hierarchy_name = getParam<std::string>("fespace_hierarchy");
83 9 : _hierarchy = getMFEMProblem().getProblemData().fespace_hierarchies.GetShared(hierarchy_name);
84 :
85 : // Parse assembly levels, optionally expanding a single input value to all levels.
86 9 : const int N = _hierarchy->GetNumLevels();
87 18 : const auto & asm_strs = getParam<std::vector<std::string>>("assembly_levels");
88 9 : const int n_asm = asm_strs.size();
89 9 : if (n_asm != 1 && n_asm != N)
90 0 : paramError(
91 : "assembly_levels", "must have length 1 or N = ", N, " (total levels), got ", n_asm, ".");
92 :
93 9 : _assembly_levels.resize(N);
94 27 : for (const auto i : make_range(N))
95 18 : _assembly_levels[i] = ParseAssemblyLevel(n_asm == 1 ? asm_strs[0] : asm_strs[i]);
96 :
97 9 : ConstructSolver();
98 9 : }
99 :
100 : void
101 9 : MFEMGeometricMultigridSolver::ConstructSolver()
102 : {
103 9 : _mg.reset();
104 9 : _level_ops.clear();
105 9 : _level_blfs.clear();
106 :
107 9 : auto proxy = std::make_unique<MGProxy>(*this);
108 9 : _mg_proxy = proxy.get();
109 9 : _solver = std::move(proxy);
110 9 : }
111 :
112 : mfem::AssemblyLevel
113 18 : MFEMGeometricMultigridSolver::ParseAssemblyLevel(const std::string & s) const
114 : {
115 18 : if (s == "legacy")
116 18 : return mfem::AssemblyLevel::LEGACY;
117 0 : if (s == "full")
118 0 : return mfem::AssemblyLevel::FULL;
119 0 : if (s == "element")
120 0 : return mfem::AssemblyLevel::ELEMENT;
121 0 : if (s == "partial")
122 0 : return mfem::AssemblyLevel::PARTIAL;
123 0 : if (s == "none")
124 0 : return mfem::AssemblyLevel::NONE;
125 0 : paramError("assembly_levels",
126 : "unknown assembly level '",
127 : s,
128 : "'. Valid values: legacy, full, element, partial, none.");
129 : return mfem::AssemblyLevel::LEGACY;
130 : }
131 :
132 : void
133 0 : MFEMGeometricMultigridSolver::SetOperator(mfem::Operator & op)
134 : {
135 0 : BuildMultigrid(op);
136 0 : }
137 :
138 : void
139 9 : MFEMGeometricMultigridSolver::BuildMultigrid(const mfem::Operator & op)
140 : {
141 9 : auto & problem = getMFEMProblem();
142 9 : auto & pd = problem.getProblemData();
143 :
144 9 : auto * eq_sys = dynamic_cast<Moose::MFEM::EquationSystem *>(pd.eqn_system.get());
145 9 : if (!eq_sys)
146 0 : mooseError("GeometricMultigridSolver '",
147 0 : name(),
148 : "': requires a standard (non-complex, non-time-dependent) EquationSystem.");
149 :
150 9 : if (eq_sys->Nonlinear())
151 2 : mooseError("GeometricMultigridSolver '",
152 2 : name(),
153 : "': nonlinear equation systems are not currently supported.");
154 :
155 7 : if (eq_sys->HasMixedBilinearForms(_var_name))
156 0 : paramError("variable",
157 : "mixed bilinear form contributions are not supported for variable '",
158 0 : _var_name,
159 : "'. Block multigrid is required for saddle-point / mixed-field problems.");
160 :
161 7 : const int N = _hierarchy->GetNumLevels();
162 7 : if (N < 1)
163 0 : paramError("fespace_hierarchy", "hierarchy must contain at least one level.");
164 7 : const int finest_level = N - 1;
165 :
166 : // Validate smoother vector length (levels 1 to N-1 each need a smoother).
167 7 : const int n_smooth = _smoother_names.size();
168 7 : if (n_smooth != 1 && n_smooth != N - 1)
169 0 : paramError("smoothers", "must have length 1 or N-1 = ", N - 1, ", got ", n_smooth, ".");
170 :
171 14 : auto get_smoother = [&](int level) -> Moose::MFEM::LinearSolverBase &
172 : {
173 14 : if (level == 0)
174 28 : return problem.getMFEMObject<Moose::MFEM::LinearSolverBase>("Moose::MFEM::SolverBase",
175 7 : _coarse_solver_name);
176 7 : const std::string & sname = (n_smooth == 1) ? _smoother_names[0] : _smoother_names[level - 1];
177 21 : return problem.getMFEMObject<Moose::MFEM::LinearSolverBase>("Moose::MFEM::SolverBase", sname);
178 7 : };
179 :
180 : // Obtain essential boundary attribute markers from the equation system.
181 7 : mfem::Array<int> ess_bdr = eq_sys->BuildEssentialBoundaryMarkers(_var_name);
182 :
183 : auto & finest_fespace =
184 7 : static_cast<mfem::ParFiniteElementSpace &>(_hierarchy->GetFESpaceAtLevel(finest_level));
185 7 : const int finest_size = finest_fespace.GetTrueVSize();
186 7 : if (op.Height() != finest_size || op.Width() != finest_size)
187 0 : mooseError("GeometricMultigridSolver '",
188 0 : name(),
189 : "': incoming fine operator has size ",
190 0 : op.Height(),
191 : " x ",
192 0 : op.Width(),
193 : ", but the finest hierarchy space has true size ",
194 : finest_size,
195 : ".");
196 :
197 : // Build new levels' forms; accumulate before touching _mg / _level_*.
198 7 : std::vector<std::shared_ptr<mfem::ParBilinearForm>> new_blfs;
199 7 : std::vector<std::unique_ptr<mfem::OperatorHandle>> new_level_ops;
200 7 : new_level_ops.reserve(N - 1);
201 :
202 7 : auto mg = std::make_unique<mfem::GeometricMultigrid>(*_hierarchy, ess_bdr);
203 7 : auto * mg_ptr = mg.get();
204 :
205 21 : for (const auto level : make_range(N))
206 : {
207 : auto & level_fespace =
208 14 : static_cast<mfem::ParFiniteElementSpace &>(_hierarchy->GetFESpaceAtLevel(level));
209 :
210 : // Compute essential true DoFs for this level.
211 14 : mfem::Array<int> level_tdofs;
212 14 : level_fespace.GetEssentialTrueDofs(ess_bdr, level_tdofs);
213 :
214 : // Build level operator.
215 14 : mfem::Operator * level_op = nullptr;
216 14 : bool own_op = false;
217 :
218 14 : if (level == finest_level)
219 : {
220 7 : level_op = const_cast<mfem::Operator *>(&op);
221 7 : own_op = false;
222 : }
223 : else
224 : {
225 : auto blf =
226 7 : eq_sys->BuildBilinearFormForFESpace(_var_name, level_fespace, _assembly_levels[level]);
227 :
228 7 : auto level_op_handle = std::make_unique<mfem::OperatorHandle>();
229 7 : blf->FormSystemMatrix(level_tdofs, *level_op_handle);
230 7 : level_op = level_op_handle->Ptr();
231 7 : own_op = false; // owned by level_op_handle or blf
232 7 : new_level_ops.push_back(std::move(level_op_handle));
233 7 : new_blfs.push_back(std::move(blf));
234 7 : }
235 :
236 : // Configure the smoother / coarse solver with this level's operator.
237 : // Each smoother's SetOperator() owns full initialization.
238 14 : auto & level_smoother = get_smoother(level);
239 14 : level_smoother.SetOperator(*level_op);
240 :
241 14 : mg_ptr->AddLevel(level_op, &level_smoother.GetSolver(), own_op, /*ownSmoother=*/false);
242 14 : }
243 :
244 : // Atomically replace:
245 : // 1. Old MG freed, dropping raw pointers into level operators.
246 : // 2. Old operator handles freed before old forms they may wrap.
247 : // 3. Proxy updated to point at the new MG and level data.
248 7 : _mg = std::move(mg);
249 7 : _level_ops = std::move(new_level_ops);
250 7 : _level_blfs = std::move(new_blfs);
251 7 : _mg_proxy->setMG(*_mg);
252 7 : }
253 : #endif
|
