|
| | __init__ (self, name, dimensions, poly_degree, unknowns_per_cell_node, solutions_per_face_node, projections_per_face_node, min_h, max_h, assembly_matrix, assembly_matrix_scaling, mass_matrix, mass_matrix_scaling, face_from_cell_projection, face_from_cell_projection_scaling, cell_from_face_projection, cell_from_face_projection_scaling, riemann_matrix, boundary_matrix, cell_relaxation, face_relaxation, approximate_system_matrix, approximate_system_matrix_scaling, integration_points, solver_tolerance, custom_init_actionset=None, custom_plot_actionset=None, custom_init_cell="", enable_plotting=True, plot_predicate="true", solver_includes="", use_smart_pointer=True, track_solution_between_cycles=False) |
| | Construct the solver.
|
| |
| | add_to_Peano4_datamodel (self, datamodel, verbose) |
| |
| | add_use_statements (self, observer) |
| | This routine should still be called even if overwritten in child class.
|
| |
| | typename (self) |
| |
| | instance_name (self) |
| | Return the name of the object that will be created for this solver.
|
| |
| | create_readme_descriptor (self) |
| |
| | add_to_abstract_solver (self, header, implementation) |
| |
| | add_to_vertex_data (self, name, size) |
| | Add in extra attributes on the fly.
|
| |
| | add_to_cell_data (self, name, size) |
| | Add in extra attributes on the fly.
|
| |
| | add_to_face_data (self, name, size) |
| | Add in extra attributes on the fly.
|
| |
| | name (self) |
| |
| | __str__ (self) |
| |
Most efficient solver we have for this type of problem class.
You can find each of the action sets here: solvers.api.actionsets.DGSolver.InitHigherOrderDofs.InitHigherOrderDofs
The text below specifies the semantics and behaviour of this matrix-free DG solver. A tutorial running step by step through the steps how to use it can be found in the tutorials section. The solver in general tackles a PDE
\( \mathcal{L}(u) = f \)
resulting in a linear equation system
\( Au = b. \)
Throughout the discussion, we unfold the compute steps over the compute mesh.
Storage
This solver holds \( K^c (p+1)^d \) unknowns (doubles) per cell for the solution, \( K^c (p+1)^d \) per cell for the right-hand side, and the same number for the residual. Holding the residual in the cell allows us to separate some algorithm steps into different action sets. In particular, we project the face solution back into the cells after updating the cell solution (in separate sweeps). These are the primary variables, i.e. those directly resulting from the PDE formulation. The correspond to the unknowns directly from \( u \) and the right-hand side.
Further to that, we hold auxiliary variables on the faces: Our Riemann solver accepts the left and the right solution. Here, each degree of freedom has \( K^f_{f \gets c} \) unknowns. It will return a result (the flux), which has \( K^f \) unknowns. Per face node, we hence hold, in total, \( ( K^f + 2K^f_{f \gets c} ) (p+1)^{d-1} \) unknowns. That is the solution for each unknown per face node, as well as a projection from the left and right for each unknown per cell. \( K^c \) and \( K^f \) can be different, but often they are the same. The user can choose how many of the \( K^c \) to project onto the face by modifying \( K^f_{f \gets c} \).
Operators
With the auxiliary variables, in place, we can write down the DG problem as
\(
\left(
\begin{array}{cccc}
A^{cc} & & & P_{c \gets f} \\
P_{f \gets c} & id & & \\
P_{f \gets c} & & id \\
& id & id & A^{ff}
\end{array}
\right)
\left(
\begin{array}{c}
u^c \\
u^+ \\
u^- \\
q
\end{array}
\right)
=
\left(
\begin{array}{c}
M b \\
0 \\
0 \\
0
\end{array}
\right)
\)
in line with xxxx. The \( q \) is the result of the Riemann solver.
Solver
We now realise a Block-Jacobi smoother subject to two relaxation parameters \( \omega _c, \omega _f \in [0,1] \). We tackle the first line in the block system with one smoothing step of a stationary iterative scheme, solve two and three exact (as they are trivial with the identity on the diagonal) and apply Jacobi to the last row again. As the first and last block are updated in parallel by the Jacobi scheme, the resuling algorithm is globally a Jacobi relaxation, too.
- Per cell, we tackle the equation \(
A^{cc} u^c + P_{c \gets f} q = M b
\) which yields the update rule
\begin{eqnarray*}
u|_c & \gets & u|_c + \omega _c \tilde {A}^{-1}\left( Mb - A^{cc} u^c - P_{c \gets f} q \right).
\end{eqnarray*}
Here, the \( q \) is the solution of the PDE term on the faces which we assume to be fixed. At this point, we are quite generic and dump the responsibility onto the user to give us an approximation of \( A_{cc} \) that we can invert. Notably, we highlight that the user may give us \( \tilde{A} = A^{cc} \) in a DG sense, as we cannot invert \( A^{cc} \).
- Per face, we use a plain Jacobi. We assume that we have an explicit expression of the Riemann solution (simple averaging would be an example) and hence \( \hat A = id \). We consequently get
\begin{eqnarray*}
q|_f & \gets & q|_f + \omega _f \left( \mathcal{R}(u^-, u^+) - q|_f \right) \\
& = & (1-\omega _f) q|_f + \omega _f \mathcal{R}(u^-, u^+).
\end{eqnarray*}
Algorithmical blueprint
- Enter cell: Whenever a cell is entered, we create a temporary \( K^c(p+1)^d \) array (on the heap).
- Into this temporary array \( r \), we write \( Mb \). This is the mass matrix (right-hand side matrix of the class constructor, as we are not restricted to a traditional mass matrix) times the right-hand sides within the cell. So \(
r \gets Mb.
\)
- Once complete, we substract \( P^{c \gets f} q|_f \). That is, we take the solution on the face and we project it (cell from face) onto the cell, where we substract it from the previous \( r \) value. We also do the same with the cell-to-cell couplings. So: \(
r \gets Mb - Pq|_f -A^{cc}u^c.
\) We also do the same with the cell couplings
- The third step is where the MFlops come from: We take the user's approximation of the cell matrix and explicitly invert it. For a linear, invariant matrix \( A \), we could do so once, but the solver itself makes no assumption if the matrix is invariant or not. So it inverts it every step and uses BLAS routines or this. The result is multiplied with the residual as accumulated so far: \(
r \gets \tilde{A}^{-1} \left( Mb - Pq|_f -A^{cc}u^c \right).
\)
- We run over the solution within the cell and scale each entry with \( (1-\omega _c) \). If \( \omega _c =1 \), this means we set the entries to zero. In general, the entries however will only be damped. After that, we add the \( r \) values to each point. This gives us a new solution guess.
- We project the updated solution \( u|_c \) onto the faces, i.e. obtain a new \( u|_f^{\pm } \). This is realised within the action set ProjectOntoFaces.
touchFaceFirstTime(): Whenever a face is touched, we solve the Riemann problem and update the Riemann guess accordingly. Initially, we assume that the \( q \) entries on the face are all set to zero.
- We evaluate per node on the face the expression \( \mathcal{R}(u^-,u^+) \). For this, we loop over the \( (p+1)^{d-1} \) points and call the C++ solver's routine. The user can thus implement the actual Riemann solver themself.
- We take the result and update the value on the face according to \(
q|_f \gets (1-\omega _f) q|_f + \omega _f \mathcal{R}(u^-, u^+).
\)
All of this functionality can be found in the subpackage actionsets.UpdateResidual.
How does it work?
Python stuff
There is a very simple Poisson solver tutorial in the tutorials section which illustrates how to use this solver.
The user needs to supply a few crucial matrices, referred to in the maths above. These are noted and their shapes discussed in the parameters section of the documentation of this class.
Under the hood
Inside the cells we store \( K^c (p+1)^{d} \) values for the solution, and the same number for the rhs and residual. We store the residual in the cell so that we can calculate the update to the cell solution separately from the projection of face solution back into the cells, but crucially are able to update the global residual before we apply the preconditioner. On each face we store \( K^f (p+1)^{d-1} \) values for the solution, and then \( 2K^f_{f \gets c} (p+1)^{d-1} \) for the projections in a separate array. We note that this factor of two comes from the fact that we have projections from cells on both the left and right. The user is free to choose how many variables from the cell they want to project onto the face. The value K^f_{f \gets c} is exposed in dg.py as the argument --projections_per_face"</tt>.
@subsubsection autotoc_md211 Matrices we pass in to the solver class
See documentation for this class itself...
So far we have three non-trivial action sets.
@subsubsection autotoc_md212 InitHigherOrderDofs
Code for this is in actionsets.InitHigherOrderDofs. So far nothing happens, other than determining which cells/faces are
on the boundary and putting in initial guesses for each cell and face dof. We pass over to the implementation
cpp file which should be left in situ in the directory for the application.
@subsubsection autotoc_md213 UpdateFaceSolution
Code for this is in actionsets.DGSolver.UpdateFaceSolution. This action set only contains something nontrivial at
touchFaceFirstTime. We apply the Riemann matrix that is present in the Abstract Solver class. This action set goes first,
the idea being that we updated the face solution using the updated face projections from the previous sweep.
@subsubsection autotoc_md214 ProjectIntoCellAndUpdateCellSolution
Code for this is in actionsets.DGSolver.ProjectIntoCellAndUpdateCellSolution. This action set works when we touch the cell for the first time.
We project the newly updated face solution back into the cells. In particular we:
- Write the updated face solution into the residual stored in the cell. This finishes the cell update from the previous sweep.
- Update the global residual
- Update the cell solution, which is projected back onto faces in the next action set.
@subsubsection autotoc_md215 ProjectOntoFaces
Code for this is in actionsets.DGSolver.ProjectOntoFaces. This occurs strictly after ProjectIntoCellAndUpdateCellSolution, so that we project
freshly updated Cell values onto the Faces.
@subsubsection autotoc_md216 UpdateResidual
Code for this is in actionsets.DGSolver.UpdateResidual. We do something nontrivial in touchCellFirstTime, namely perform the cell update
rule as above. We set the residual to be \_form#322.
@param solver_includes C++ code snippet
Arbitrary C++ code snippet. Usually used to include additional files
(therefore the name), though you can use it for preprocessing macros
or using statements as well. Has to be a valid C++ expression, i.e.
something similar to
@icode
#include "myfile.h"
Definition at line 21 of file DiscontinuousGalerkinDiscretisationPointWiseRiemannSolver.py.
Public Member Functions inherited from solvers.api.Solver.Solver
1.9.8