Stationary Solvers

Meshes and FESpaces

To solve a ProblemDescription the user needs to provide some discretisation information. The most important one is the mesh (that should be conforming with the region information used in the problem description).

Meshes are expected in the form of an ExtendableGrid, see ExtendableGrids.jl for details and constructors. Grid generators for simplex grids can be found e.g. in the external module SimplexGridFactory.jl. There is also an extension to read meshes from gmsh files.

Secondly, the user needs to decide how to discretize the unknowns by providing a FESpace for each unknown, which can be generated by

FESpace{FEType}(grid::ExtendableGrid)

Here, FEType denotes the type of finite element. A list of available FETypes can be found in the Documentation of ExtendableFEMBase.jl.

Solve (monolithic)

If solve is applied to a PDEDescription and an array of FESpaces (that specify the ansatz spaces for the unknowns) an investigation of the PDEDescription is performed that decides if the problem is nonlinear (and has to be solved by a fixed-point algorithm) or if it can be solved directly in one step. The solver attempts to bring the (nonlinear) residual to the prescribed target tolerance.

CommonSolve.solve — Function
function solve(
	PD::ProblemDescription,
	[FES::Union{<:FESpace,Vector{<:FESpace}}],
	SC = nothing;
	unknowns = PD.unknowns,
	kwargs...)

Returns a solution of the PDE as an FEVector for the provided FESpace(s) FES (to be used to discretised the unknowns of the PDEs). If no FES is provided an initial FEVector (see keyword init) must be provided (which is used to built the FES).

This function extends the CommonSolve.solve interface and the PDEDescription takes the role of the ProblemType and FES takes the role of the SolverType.

Keyword arguments:

  • verbosity: verbosity level. Default: 0

  • method_linear: any solver or custom LinearSolveFunction compatible with LinearSolve.jl (default = UMFPACKFactorization()). Default: LinearSolve.UMFPACKFactorization(true, true)

  • time: current time to be used in all time-dependent operators. Default: 0.0

  • spy: show unicode spy plot of system matrix during solve. Default: false

  • show_config: show configuration at the beginning of solve. Default: false

  • preconlinear: function that computes preconditioner for methodlinear incase an iterative solver is chosen. Default: nothing

  • init: initial solution (also used to save the new solution). Default: nothing

  • initialized: linear system in solver configuration is already assembled (turns true after first solve). Default: false

  • return_config: solver returns solver configuration (including A and b of last iteration). Default: false

  • show_matrix: show system matrix after assembly. Default: false

  • reltol: reltol for linear solver (if iterative). Default: 1.0e-11

  • is_linear: linear problem (avoid reassembly of nonlinear operators to check residual). Default: ''auto''

  • symmetrize: make system matrix symmetric (replace by (A+A')/2). Default: false

  • inactive: inactive unknowns (are made available in assembly, but not updated in solve). Default: ExtendableFEM.Unknown[]

  • plot: plot all solved unknowns with a (very rough but fast) unicode plot. Default: false

  • damping: amount of damping, value should be between in (0,1). Default: 0

  • check_matrix: check matrix for symmetry and positive definiteness and largest/smallest eigenvalues. Default: false

  • restrict_dofs: array of dofs for each unknown that should be solved (default: all dofs). Default: Any[]

  • constant_rhs: right-hand side is constant (skips reassembly). Default: false

  • abstol: abstol for linear solver (if iterative). Default: 1.0e-11

  • symmetrize_structure: make the system sparse matrix structurally symmetric (e.g. if [j,k] is also [k,j] must be set, all diagonal entries must be set). Default: false

  • target_residual: stop if the absolute (nonlinear) residual is smaller than this number. Default: 1.0e-10

  • maxiterations: maximal number of nonlinear iterations/linear solves. Default: 10

  • constant_matrix: matrix is constant (skips reassembly and refactorization in solver). Default: false

Depending on the detected/configured nonlinearities the whole system is either solved directly in one step or via a fixed-point iteration.

source
Note

The type of fixed-point algorithm depends on the discretisation of the nonlinearities. If all of them are assembled as a NonlinearOperator, this will result in a Newton scheme (which can be somewhat costumized via the keywords arguments like damping). If all nonlinearities are linearized by LinearOperator and BilinearOperator, this will result in other types of fixed-point iterations.

Solve (iterating subproblems)

If the problem can be solved by iterating over subproblems this can be achieved as well. For that each subproblem must be configured separately via a SolverConfiguration (allowing different tolerances and keyword arguments for each subproblem solve). A SolverConfiguration can be constructed with this constructor:

ExtendableFEM.SolverConfiguration — Type
function iterate_until_stationarity(
	SolverConfiguration(Problem::ProblemDescription
	[FES::Union{<:FESpace, Vector{<:FESpace}}];
	init = nothing,
	unknowns = Problem.unknowns,
	kwargs...)

Returns a solver configuration for the ProblemDescription that can be passed to the solve function. Here, FES are the FESpaces that should be used to discretize the selected unknowns. If no FES is provided an initial FEVector (see keyword init) must be provided (which is used to built the FES).

Keyword arguments:

  • verbosity: verbosity level. Default: 0

  • method_linear: any solver or custom LinearSolveFunction compatible with LinearSolve.jl (default = UMFPACKFactorization()). Default: LinearSolve.UMFPACKFactorization(true, true)

  • time: current time to be used in all time-dependent operators. Default: 0.0

  • spy: show unicode spy plot of system matrix during solve. Default: false

  • show_config: show configuration at the beginning of solve. Default: false

  • preconlinear: function that computes preconditioner for methodlinear incase an iterative solver is chosen. Default: nothing

  • init: initial solution (also used to save the new solution). Default: nothing

  • initialized: linear system in solver configuration is already assembled (turns true after first solve). Default: false

  • return_config: solver returns solver configuration (including A and b of last iteration). Default: false

  • show_matrix: show system matrix after assembly. Default: false

  • reltol: reltol for linear solver (if iterative). Default: 1.0e-11

  • is_linear: linear problem (avoid reassembly of nonlinear operators to check residual). Default: ''auto''

  • symmetrize: make system matrix symmetric (replace by (A+A')/2). Default: false

  • inactive: inactive unknowns (are made available in assembly, but not updated in solve). Default: ExtendableFEM.Unknown[]

  • plot: plot all solved unknowns with a (very rough but fast) unicode plot. Default: false

  • damping: amount of damping, value should be between in (0,1). Default: 0

  • check_matrix: check matrix for symmetry and positive definiteness and largest/smallest eigenvalues. Default: false

  • restrict_dofs: array of dofs for each unknown that should be solved (default: all dofs). Default: Any[]

  • constant_rhs: right-hand side is constant (skips reassembly). Default: false

  • abstol: abstol for linear solver (if iterative). Default: 1.0e-11

  • symmetrize_structure: make the system sparse matrix structurally symmetric (e.g. if [j,k] is also [k,j] must be set, all diagonal entries must be set). Default: false

  • target_residual: stop if the absolute (nonlinear) residual is smaller than this number. Default: 1.0e-10

  • maxiterations: maximal number of nonlinear iterations/linear solves. Default: 10

  • constant_matrix: matrix is constant (skips reassembly and refactorization in solver). Default: false

source

If each subproblem has a SolverConfiguration the fixed-point iteration can be triggered with this function:

ExtendableFEM.iterate_until_stationarity — Function
function iterate_until_stationarity(
	SCs::Array{<:SolverConfiguration, 1},
	FES = nothing;
	maxsteps = 1000,
	init = nothing,
	unknowns = [SC.PD.unknowns for SC in SCs],
	kwargs...)

Iterates consecutively over all SolverConfigurations (each contains the ProblemDescription of the corressponding subproblem) until the residuals of all subproblems are below their tolerances and returns the solution of the combined unknowns of all subproblems. The additional argument maxsteps limits the number of these iterations If an initial vector init is provided it should contain all unknowns of the subproblems.

Using the SolverConfiguration instead of the ProblemDescription in the first argument allows to use different kwargs for each subproblem. The SolverConfiguration for each subproblem can be generated by

SolverConfiguration(PD::ProblemDescription; init = sol, kwargs...)

with the usual keyword arguments.

source