20 #include "libmesh/libmesh_common.h" 21 #include "libmesh/equation_systems.h" 22 #include "libmesh/nonlinear_implicit_system.h" 23 #include "libmesh/nonlinear_solver.h" 24 #include "libmesh/linear_implicit_system.h" 25 #include "libmesh/transient_system.h" 26 #include "libmesh/numeric_vector.h" 27 #include "libmesh/sparse_matrix.h" 28 #include "libmesh/string_to_enum.h" 29 #include "libmesh/coupling_matrix.h" 38 params.
addClassDescription(
"Physics-based preconditioner (PBP) allows individual physics to have " 39 "their own preconditioner.");
43 "The order the block rows will be solved in. Put the name of variables here " 44 "to stand for solving that variable's block row. A variable may appear more " 45 "than once (to create cylces if you like).");
46 params.
addRequiredParam<std::vector<std::string>>(
"preconditioner",
"TODO: docstring");
54 _nl(_fe_problem.getNonlinearSystemBase(_nl_sys_num))
56 const auto & libmesh_system =
_nl.
system();
57 unsigned int num_systems =
_nl.
system().n_vars();
69 std::unique_ptr<CouplingMatrix> cm = std::make_unique<CouplingMatrix>(
n_vars);
71 bool full = getParam<bool>(
"full");
76 for (
unsigned int i = 0; i <
n_vars; i++)
80 for (
const auto & off_diag : getParam<NonlinearVariableName, NonlinearVariableName>(
81 "off_diag_row",
"off_diag_column"))
83 const auto row = libmesh_system.variable_number(off_diag.first);
84 const auto column = libmesh_system.variable_number(off_diag.second);
85 (*cm)(row, column) = 1;
92 for (
unsigned int i = 0; i <
n_vars; i++)
93 for (
unsigned int j = 0; j <
n_vars; j++)
101 const std::vector<std::string> & pc_types = getParam<std::vector<std::string>>(
"preconditioner");
102 for (
unsigned int i = 0; i < num_systems; i++)
103 _pre_type[i] = Utility::string_to_enum<PreconditionerType>(pc_types[i]);
106 const std::vector<std::string> & solve_order = getParam<std::vector<std::string>>(
"solve_order");
115 std::vector<std::vector<unsigned int>> off_diag(
n_vars);
118 const std::vector<NonlinearVariableName> & odr =
119 getParam<std::vector<NonlinearVariableName>>(
"off_diag_row");
120 const std::vector<NonlinearVariableName> & odc =
121 getParam<std::vector<NonlinearVariableName>>(
"off_diag_column");
125 unsigned int row =
_nl.
system().variable_number(odr[i]);
126 unsigned int column =
_nl.
system().variable_number(odc[i]);
127 off_diag[row].push_back(column);
131 for (
unsigned int var = 0; var <
n_vars; var++)
137 mooseError(
"PBP must be used with JFNK solve type");
144 std::vector<unsigned int> off_diag,
147 std::string var_name =
_nl.
system().variable_name(var);
149 LinearImplicitSystem & precond_system =
150 _fe_problem.
es().add_system<LinearImplicitSystem>(var_name +
"_system");
151 precond_system.assemble_before_solve =
false;
154 precond_system.add_variable(
155 var_name +
"_prec",
_nl.
system().variable(var).type(), active_subdomains);
173 TIME_SECTION(
"init", 2,
"Initializing PhysicsBasedPreconditioner");
178 const unsigned int num_systems =
_systems.size();
184 for (
unsigned int i = 0; i < num_systems; i++)
189 for (
unsigned int system_var = 0; system_var < num_systems; system_var++)
191 LinearImplicitSystem & u_system = *
_systems[system_var];
200 preconditioner->
set_matrix(u_system.get_system_matrix());
203 preconditioner->
init();
210 const unsigned int num_systems =
_systems.size();
212 std::vector<JacobianBlock *> blocks;
215 for (
unsigned int system_var = 0; system_var < num_systems; system_var++)
217 LinearImplicitSystem & u_system = *
_systems[system_var];
221 new JacobianBlock(u_system, u_system.get_system_matrix(), system_var, system_var);
222 blocks.push_back(block);
227 unsigned int coupled_var =
_off_diag[system_var][diag];
228 std::string coupled_name =
_nl.
system().variable_name(coupled_var);
232 blocks.push_back(block);
239 for (
auto & block : blocks)
246 TIME_SECTION(
"apply", 1,
"Applying PhysicsBasedPreconditioner");
248 const unsigned int num_systems =
_systems.size();
253 for (
unsigned int sys = 0; sys < num_systems; sys++)
261 LinearImplicitSystem & u_system = *
_systems[system_var];
265 mesh,
_nl.
system().number(), system_var, x, u_system.number(), 0, *u_system.rhs);
271 unsigned int coupled_var =
_off_diag[system_var][diag];
272 LinearImplicitSystem & coupled_system = *
_systems[coupled_var];
290 u_system.rhs->close();
300 for (
unsigned int system_var = 0; system_var < num_systems; system_var++)
302 LinearImplicitSystem & u_system = *
_systems[system_var];
305 mesh, u_system.number(), 0, *u_system.solution,
_nl.
system().number(), system_var, y);
virtual void setup()
This is called every time the "operator might have changed".
virtual const std::set< SubdomainID > * getVariableBlocks(unsigned int var_number)
Get the block where a variable of this system is defined.
Helper class for holding the preconditioning blocks to fill.
virtual void apply(const NumericVector< Number > &x, NumericVector< Number > &y)
Computes the preconditioned vector "y" based on input "x".
static void copyVarValues(MeshBase &mesh, const unsigned int from_system, const unsigned int from_var, const NumericVector< Number > &from_vector, const unsigned int to_system, const unsigned int to_var, NumericVector< Number > &to_vector)
Helper function for copying values associated with variables in vectors from two different systems...
std::vector< std::vector< unsigned int > > _off_diag
Holds which off diagonal blocks to compute.
SolverParams & solverParams()
Get the solver parameters.
virtual void attachPreconditioner(Preconditioner< Number > *preconditioner)=0
Attach a customized preconditioner that requires physics knowledge.
PreconditionerType type() const
NonlinearSystemBase & _nl
The nonlinear system this PBP is associated with (convenience reference)
std::vector< PreconditionerType > _pre_type
Which preconditioner to use for each solve.
std::vector< LinearImplicitSystem * > _systems
List of linear system that build up the preconditioner.
void set_matrix(SparseMatrix< Number > &mat)
FEProblemBase & _fe_problem
Subproblem this preconditioner is part of.
const Parallel::Communicator & _communicator
virtual void clear()
Release all memory and clear data structures.
Base class for MOOSE preconditioners.
virtual EquationSystems & es() override
virtual void init()
Initialize data structures if not done so already.
bool isParamValid(const std::string &name) const
Test if the supplied parameter is valid.
virtual unsigned int nVariables() const
Get the number of variables in this system.
virtual void computeJacobianBlocks(std::vector< JacobianBlock *> &blocks, const unsigned int nl_sys_num)
Computes several Jacobian blocks simultaneously, summing their contributions into smaller preconditio...
virtual void scale(const Number factor)=0
static InputParameters validParams()
Constructor.
virtual ~PhysicsBasedPreconditioner()
PhysicsBasedPreconditioner(const InputParameters ¶ms)
Jacobian-Free Newton Krylov.
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Implements a segregated solve preconditioner.
std::vector< std::vector< SparseMatrix< Number > * > > _off_diag_mats
Holds pointers to the off-diagonal matrices.
void addSystem(unsigned int var, std::vector< unsigned int > off_diag, PreconditionerType type=AMG_PRECOND)
Add a diagonal system + possibly off-diagonals ones as well, also specifying type of preconditioning...
static InputParameters validParams()
unsigned int number() const
Gets the number of this system.
void vector_mult_add(NumericVector< Number > &dest, const NumericVector< Number > &arg) const
virtual System & system() override
Get the reference to the libMesh system.
std::vector< std::unique_ptr< Preconditioner< Number > > > _preconditioners
Holds one Preconditioner object per small system to solve.
virtual MooseMesh & mesh() override
std::vector< unsigned int > _solve_order
Holds the order the blocks are solved for.
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
void set_type(const PreconditionerType pct)
static std::unique_ptr< Preconditioner< Number > > build_preconditioner(const libMesh::Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
registerMooseObjectAliased("MooseApp", PhysicsBasedPreconditioner, "PBP")
void setCouplingMatrix(std::unique_ptr< CouplingMatrix > cm)
Setup the coupling matrix on the finite element problem.
auto index_range(const T &sizable)