10from .Solver
import Solver
11from multigrid.petsc.api.actionsets.EnumerateDoFs
import EnumerateDoFs
12from multigrid.petsc.api.actionsets.InitCellDoFs
import InitCellDoFs
13from multigrid.petsc.api.actionsets.ProjectPETScSolutionBackOntoMesh
import ProjectPETScSolutionOnCellsBackOntoMesh
14from multigrid.petsc.api.actionsets.PlotDGDataInPeanoBlockFormat
import PlotDGDataInPeanoBlockFormat
15from multigrid.petsc.api.actionsets.PlotExactSolution
import PlotExactSolution
16from multigrid.petsc.api.actionsets.ImposeDirichletBoundaryConditionsWithInteriorPenaltyMethod
import ImposeDirichletBoundaryConditionsWithInteriorPenaltyMethod
18from abc
import abstractmethod
27 Trigger assembly of PETSc matrix and also assemble rhs vector
30 TemplateTouchCellFirstTime =
"""
31 if (fineGridCell{{SOLVER_NAME}}PETScData.getType() == celldata::{{SOLVER_NAME}}PETScData::Type::Interior) {
32 std::pair<int,int> localCellIndex = std::make_pair(_spacetreeId, fineGridCell{{SOLVER_NAME}}PETScData.getUnknownBaseNumber());
33 int globalCellIndex = repositories::{{SOLVER_INSTANCE}}.getLocalToGlobalMap().getGlobalIndex(localCellIndex, ::petsc::LocalToGlobalMap::Type::Cell);
35 // ================================
36 // Left-hand side A_cc contribution
37 // ================================
38 auto lhsCellMatrix = repositories::{{SOLVER_INSTANCE}}.getLhsMatrix(marker.x(),marker.h());
40 logTraceInWith3Arguments( "touchCellFirstTime()", globalCellIndex, lhsCellMatrix.toString(), marker.toString() );
42 for (int row=0; row<repositories::{{SOLVER_INSTANCE}}.DoFsPerCell; row++)
43 for (int col=0; col<repositories::{{SOLVER_INSTANCE}}.DoFsPerCell; col++) {
44 repositories::{{SOLVER_INSTANCE}}.getLinearEquationSystem().insert(
47 lhsCellMatrix(row,col)
50 logTraceOut( "touchCellFirstTime()" );
52 // ================================
53 // Right-hand side vector contribution (multiply rhs (mass) matrix with rhs initial values
54 // ================================
55 auto rhsCellMatrix = repositories::{{SOLVER_INSTANCE}}.getRhsMatrix(marker.x(),marker.h());
57 for (int row=0; row<repositories::{{SOLVER_INSTANCE}}.DoFsPerCell; row++) {
59 for (int col=0; col<repositories::{{SOLVER_INSTANCE}}.DoFsPerCell; col++) {
60 rhs += rhsCellMatrix(row,col) * fineGridCell{{SOLVER_NAME}}.getRhs(col);
62 repositories::{{SOLVER_INSTANCE}}.getLinearEquationSystem().insert(
63 globalCellIndex+row, rhs
67 auto projectionCellToFaceMatrix = repositories::{{SOLVER_INSTANCE}}.getProjectionOfCellDataOntoFace(marker.x(),marker.h());
68 auto projectionFaceToCellMatrix = repositories::{{SOLVER_INSTANCE}}.getProjectionOfRiemannSolutionOntoCell(marker.x(),marker.h());
69 int dofsPerFace = repositories::{{SOLVER_INSTANCE}}.NodesPerFace*repositories::{{SOLVER_INSTANCE}}.FaceUnknowns;
71 for (int d=0; d<Dimensions; d++) {
72 // =================================
73 // Project to left face along axis d
74 // =================================
76 fineGridFaces{{SOLVER_NAME}}PETScData(d).getType()==facedata::{{SOLVER_NAME}}PETScData::Type::Interior
78 fineGridFaces{{SOLVER_NAME}}PETScData(d).getType()==facedata::{{SOLVER_NAME}}PETScData::Type::Boundary
80 std::pair<int,int> localFirstDoFIndexOfLeftFace = std::make_pair(_spacetreeId, fineGridFaces{{SOLVER_NAME}}PETScData(d).getUnknownBaseNumber());
81 int globalFirstDoFIndexOfLeftFace = repositories::{{SOLVER_INSTANCE}}.getLocalToGlobalMap().getGlobalIndex(localFirstDoFIndexOfLeftFace, ::petsc::LocalToGlobalMap::Type::Face);
83 assertion1(globalFirstDoFIndexOfLeftFace>=0, marker.toString());
84 logTraceInWith5Arguments( "touchCellFirstTime::LeftProjection", globalFirstDoFIndexOfLeftFace, globalCellIndex , projectionCellToFaceMatrix, projectionFaceToCellMatrix, marker.toString() );
86 for (int row=0; row<dofsPerFace; row++ ) {
87 for (int col=0; col<repositories::{{SOLVER_INSTANCE}}.DoFsPerCell; col++) {
88 logTraceIn("CellToFace");
89 repositories::{{SOLVER_INSTANCE}}.getLinearEquationSystem().insert(
90 globalFirstDoFIndexOfLeftFace + row + dofsPerFace,
92 projectionCellToFaceMatrix(d * dofsPerFace + row,col)
94 logTraceOut("CellToFace");
96 logTraceInWith1Argument("FaceToCellInsertion", projectionFaceToCellMatrix(col, d * dofsPerFace + row));
97 repositories::{{SOLVER_INSTANCE}}.getLinearEquationSystem().insert(
99 globalFirstDoFIndexOfLeftFace + row + 2 * dofsPerFace,
100 projectionFaceToCellMatrix(col, d * dofsPerFace + row)
102 logTraceOut("FaceToCell");
105 logTraceOut( "touchCellFirstTime::LeftProjection" );
109 fineGridFaces{{SOLVER_NAME}}PETScData(d+Dimensions).getType()==facedata::{{SOLVER_NAME}}PETScData::Type::Interior
111 fineGridFaces{{SOLVER_NAME}}PETScData(d+Dimensions).getType()==facedata::{{SOLVER_NAME}}PETScData::Type::Boundary
113 std::pair<int,int> localFirstDoFIndexOfRightFace = std::make_pair(_spacetreeId, fineGridFaces{{SOLVER_NAME}}PETScData(d+Dimensions).getUnknownBaseNumber());
114 int globalFirstDoFIndexOfRightFace = repositories::{{SOLVER_INSTANCE}}.getLocalToGlobalMap().getGlobalIndex(localFirstDoFIndexOfRightFace, ::petsc::LocalToGlobalMap::Type::Face);
116 assertion1(globalFirstDoFIndexOfRightFace>=0, marker.toString());
118 logTraceInWith2Arguments( "touchCellFirstTime::RightProjection", globalFirstDoFIndexOfRightFace, marker.toString() );
120 for (int row=0; row<dofsPerFace; row++ ) {
121 for (int col=0; col<repositories::{{SOLVER_INSTANCE}}.DoFsPerCell; col++) {
122 repositories::{{SOLVER_INSTANCE}}.getLinearEquationSystem().insert(
123 globalFirstDoFIndexOfRightFace + row,
125 projectionCellToFaceMatrix(d * dofsPerFace + row + Dimensions*dofsPerFace,col)
127 repositories::{{SOLVER_INSTANCE}}.getLinearEquationSystem().insert(
129 globalFirstDoFIndexOfRightFace + row + 2 * dofsPerFace,
130 projectionFaceToCellMatrix(col, d * dofsPerFace + row + Dimensions*dofsPerFace)
134 } // end of if interior face
136 logTraceOut( "touchCellFirstTime::RightProjection" );
142 TemplateTouchFaceFirstTime=
"""
143 auto RiemannSolverMatrix = repositories::{{SOLVER_INSTANCE}}.getRiemannSolver(marker.x(),marker.h());
144 int dofsPerFace = repositories::{{SOLVER_INSTANCE}}.NodesPerFace*repositories::{{SOLVER_INSTANCE}}.FaceUnknowns;
145 if ( fineGridFace{{SOLVER_NAME}}PETScData.getType() == facedata::{{SOLVER_NAME}}PETScData::Type::Interior ) {
146 std::pair<int,int> localFirstDoFIndex = std::make_pair(_spacetreeId, fineGridFace{{SOLVER_NAME}}PETScData.getUnknownBaseNumber());
147 int globalFirstDoFIndex = repositories::{{SOLVER_INSTANCE}}.getLocalToGlobalMap().getGlobalIndex(localFirstDoFIndex, ::petsc::LocalToGlobalMap::Type::Face);
149 assertion( globalFirstDoFIndex>=0 );
151 logTraceInWith2Arguments("touchFaceFirstTime", globalFirstDoFIndex, RiemannSolverMatrix);
152 for (int row=0; row<dofsPerFace; row++)
154 // we have col<2*dofsPerFace here since we have left and right projections here.
155 // ie, we have a certain number of "true" dofsPerFace and then twice as many
156 // helper dofs to hold the projections.
157 // In this loop, we are setting the values of the true face Dofs, using the
158 // left and right projections. This is why we get the global index of the
159 // face and offset by 2*dofsPerFace: this will index the true dofs,
160 // which are enumerated AFTER the projections.
161 for (int col=0; col<2*dofsPerFace; col++) {
162 repositories::{{SOLVER_INSTANCE}}.getLinearEquationSystem().insert(
163 globalFirstDoFIndex + 2*dofsPerFace + row,
164 globalFirstDoFIndex + col,
165 RiemannSolverMatrix(row, col)
169 // 3*dofsPerFace so that each of the -, + and solution variables receive an identity here.
170 for (int i=0; i<3*dofsPerFace; i++)
172 repositories::{{SOLVER_INSTANCE}}.getLinearEquationSystem().insert(
173 globalFirstDoFIndex + i,
174 globalFirstDoFIndex + i,
179 logTraceOut("touchFaceFirstTime");
182 else if ( fineGridFace{{SOLVER_NAME}}PETScData.getType() == facedata::{{SOLVER_NAME}}PETScData::Type::Boundary ){
183 // add some more identities!
196 super( AssemblePetscMatrix, self ).
__init__()
198 self.
d[
"SOLVER_INSTANCE"] = solver.instance_name()
199 self.
d[
"SOLVER_NAME"] = solver.typename()
205Provide C++ code snippet for peano4.solversteps.ActionSet.OPERATION_TOUCH_VERTEX_FIRST_TIME
206Provide C++ code snippet for peano4.solversteps.ActionSet.OPERATION_TOUCH_CELL_FIRST_TIME
208Only touchVertexFirstTime is an event where this action set actually
209does something: It inserts the template TemplateTouchVertexFirstTime and
210replaces it with entries from the dictionary. The latter is befilled
213We actually do something during touchVertexFirstTime and touchCellFirstTime. We insert the
214appropriate template into each.
218 if operation_name==peano4.solversteps.ActionSet.OPERATION_TOUCH_CELL_FIRST_TIME:
221 if operation_name==peano4.solversteps.ActionSet.OPERATION_TOUCH_FACE_FIRST_TIME:
230 Configure name of generated C++ action set
232 This action set will end up in the directory observers with a name that
233 reflects how the observer (initialisation) is mapped onto this action
234 set. The name pattern is ObserverName2ActionSetIdentifier where this
235 routine co-determines the ActionSetIdentifier. We make is reflect the
239 return __name__.replace(
".py",
"").replace(
".",
"_")
245 The action set that Peano will generate that corresponds to this class
246 should not be modified by users and can safely be overwritten every time
247 we run the Python toolkit.
256Consult petsc.Project for details
260#include "repositories/SolverRepository.h"
261#include "tarch/la/Matrix.h"
276Define body of constructor
282 _spacetreeId = treeNumber;
287class DiscontinuousGalerkinDiscretisation(
Solver):
290 \page petsc_dg_solver Discontinuous Galerkin Solver with PETSc
292 ## The weak formulation
294 To translate a PDE into an equation system, we express the solution
296 @f$ u = \sum _i u_i \phi _i @f$
298 as sum over shape functions @f$ \phi _i @f$.
299 The shape functions in this class are locally continuous, but their linear
300 combination is not. We test this setup against test functions @f$ v @f$.
301 Each test function is actually any of the shape functions @f$ \phi _i @f$.
302 Drawing both shapes and test functions from the same function space means that
303 we follow a ***Ritz-Galerkin*** approach.
306 Let @f$ \mathcal{L}(u) @f$ be the differential operator on the left-hand side
307 of the PDE, i.e. let's assume that we solve
313 A weak formulation means that we multiple both sides with a test function and
314 then integrate over the whole thing. As we take all these tests from the same
315 function space as the shape functions, we know that they all have local
316 support.Consequently, our test integral "shrinks" down to the support of the
320 \int _\Omega \mathcal{L}(u) \cdot v dx = \int_{\text{cell}_v} \mathcal{L}(u) \cdot v dx
323 as the product disappears everywhere outside of the cell that hosts the test
324 function of interest.
327 Once again, we can integrate by parts to get rid of the highest derivative in
328 our partial differential operator @f$ \mathcal{L}(u) @f$, but this time, we will
329 get an expression of the form
333 \int _\Omega \mathcal{L}(u) \cdot v dx
335 - \int _\Omega \left( \tilde{\mathcal{L}}(u), \nabla v \right) dx
336 + \int _{\partial \text{cell}_v} (n,\tilde{\mathcal{M}}(u)) v dS(x)
339 - \int _{\text{cell}_v} \left( \tilde{\mathcal{L}}(u), \nabla v \right) dx
340 + \int _{\partial \text{cell}_v} (n,\tilde{\mathcal{M}}(\hat u )) v dS(x)
342 & = & \int _{\text{cell}_v} f \cdot v dx.
345 We reduced the order of the deriative within @f$ \mathcal{L}(u) @f$ and got a
346 new operator @f$ \tilde{\mathcal{L}} @f$.
347 If @f$ \mathcal{L}(u) @f$ is dominated by a second derivative, @f$ \tilde{\mathcal{L}} @f$
348 hosts only first-order derivatives.
349 But this is a Pyrrhic victory:
350 We now get an additional PDE term over the cell faces.
355 1. the jump terms along the faces do not disappear, i.e. they do not cancel out
356 from left and right as they do in a Continuous Galerkin formulation.
357 2. Even worse, they refer to a value @f$ u @f$ which technically does not
358 exist, as there is no (unique) solution along the face. We only have a
359 left @f$ u^- @f$ and a right @f$ u^+ @f$ solution. Therefore, I explicitly
360 write @f$ \hat u @f$ in the equation above. It denotes that we haven't yet
361 decided what the function looks like along the face.
362 3. We now have cell integrals to evaluate, and there are additional face
363 integrals. The latter ones couple neighbouring cells.
365 We end up with a linear equation system
369 where x has @f$ |\mathbb{C}| \cdot K \cdot (p+1)^d @f$ entries. @f$ \mathbb{C} @f$
370 is set of cells in the mesh.
373 ## An operator language and temporary variables
375 Before we dive into examples or discuss how to map this formalism onto code, we
376 introduce an operator language. That is, we break down the overall computation
377 into small tasks. These tasks will be found 1:1 in the code later on. Furthermore,
378 we gain some freedom by breaking things down into tiny bits and pieces: We can
379 alter individual ones or puzzle them together in various ways.
381 We introduce a few operators (functions) to break down the overall computational
382 scheme into a sequence/set of computations:
384 - The operator @f$ A_{cc} @f$ takes the dofs within a cell and computes the
385 term @f$ - \int _{\text{cell}_v} \left( \tilde{\mathcal{L}}(u), \nabla v \right) dx @f$.
386 It takes the @f$ K (p+1)^d @f$ values per cell and yields @f$ K (p+1)^d @f$
387 equations/updates from the cell. In classes like petsc.solvers.DiscontinuousGalerkinDiscretisation
388 this matrix is called cell_cell_lhs_matrix.
389 - Let's introduce temporary variables on the face. Those are the
390 @f$ u^+ @f$ and @f$ u^- @f$ values, i.e. we store @f$ 2K (p+1)^{d-1} @f$
391 degrees of freedom on the face. Those are not really degrees of freedom.
392 They will be mere projection, i.e. simply encode what the solution along
394 - The operator @f$ A_{cf} @f$ takes the solution from a cell and projects it
395 onto a face. Let @f$ u^+ @f$ and @f$ u^- @f$ denote the solution left or
396 right of a face (from the face's point of view). For each of the @f$ 2d @f$
397 faces of a cell, the operator @f$ P_{cf} @f$ yields either the
398 @f$ u^+ @f$ or @f$ u^- @f$.
399 - The operator @f$ A_{ff} @f$ operates on a face. It takes
400 @f$ u^+ @f$ and @f$ u^- @f$ and yields something
401 that "reconstructs" the solution and applies the flux operator on it.
402 This is the Riemann solver and hence linear algebraises the
403 @f$ \tilde{\mathcal{M}}(\hat u ) @f$ from above.
404 - The operator @f$ A_{fc} @f$ takes the @f$ \hat u @f$ on the face and yields
405 the updates for the cell. It evaluates @f$ \int _{\partial \text{cell}_v} (n,\tilde{\mathcal{M}}(\hat u )) v dS(x) @f$.
408 <div style="background-color: #cfc ; padding: 10px; border: 1px solid green;">
409 Within Peano, it is important to distinguish data on cell and faces carefully
410 from each other. Furthermore
411 (i) we can read data associated with faces within cells;
412 (ii) we can never read data from neighbouring cells within a cell;
413 (iii) we can never access adjacent cell data from a face.
418 The overall scheme in operator notation thus reads
421 A_{cc} + A_{fc} \circ A_{ff} \circ A_{cf}
425 i.e. we evaluate the cell contributions and add the face contributions
426 which in turn results from a concatenation of three operators:
427 Project solution onto face, solve a Riemann problem there, and then
428 take the result and project it back onto the cell. The projection onto
429 the face is a purely geoemtric thing: Once we know the polynomials
430 chosen within the cell, we can evaluate what the solution looks like
431 depending on the weights (unknowns within the degree of freedom).
432 With that knowledge, we can determine what the solution on the face
433 looks like. We use the same polynomial space there (restricted to the
434 submanifold) and therefore can express the face representation through
435 polynomial weights (left and right of the face), too. For each pair of
436 the polynomial weights, we can solve the Riemann problem and store the
437 outcome in further weights. That is, we represent the outcome of the
438 Riemann problem again in the polynomial space on the faces. Finally,
439 we take this representation, test it against the test functions,
440 integrate over the face and write the contributions back into the cell.
445 Effectively, our approach blows up the equation system to solve into
450 A_{cc} & A_{fc} & 0 & 0 \\
451 0 & id & -A_{ff}|_+ & -A_{ff}|_- \\
452 -id & 0 & A_{cf}|_+ & 0 \\
453 -id & 0 & 0 & A_{cf}|_-
456 x_{c} \\ x_f \\ x^+ \\ x^-
465 @todo We need a nice illustration here of all the function spaces
467 where we have to specify the involved matrices further.
468 The @f$ x_c @f$ are the real degrees of freedom, i.e. the weights within
469 the cells which effectively span the solution.
470 The @f$ x_\pm @f$ encode the projection of these polynomials onto the faces.
471 These are not real degrees of freedom, as they always result directly from
473 The @f$ x_f @f$ encode the reconstructed solution on the faces and hold the
474 average of the @f$ x_\pm @f$ values.
475 This setup means our total matrix is no longer square.
476 @f$ A_{cc} @f$ is square, with @f$ |\mathbb{C}| \cdot (p+1)^d @f$ rows,
477 where @f$ |\mathbb{C}| @f$ is the number of cells in the mesh.
478 @f$ A_{fc} @f$ has the same number of rows, but @f$ |\mathbb{F}| \cdot (p+1)^{d-1} @f$
479 columns, since this part has to "hit" @f$ x_f @f$.
480 @f$ A_{ff} @f$ has @f$ |\mathbb{F}| \cdot (p+1)^{d-1} @f$ rows, and twice as many
481 columns (note that @f$ | x_\pm| \ = \ 2 |x_f| @f$ ).
482 Finally, @f$ A_{cf} @f$ has @f$ |\mathbb{C}| \cdot (p+1)^d @f$ rows and
483 @f$ 2|\mathbb{F}| \cdot (p+1)^{d-1} @f$ columns, since this part "hits"
486 Our global matrix is no longer square. We have
487 @f$ 2 |\mathbb{C}| \cdot (p+1)^d + |\mathbb{F}| \cdot (p+1)^{d-1} @f$ rows
489 @f$ |\mathbb{C}| \cdot (p+1)^d + 3|\mathbb{F}| \cdot (p+1)^{d-1} @f$ columns.
494 <div style="background-color: #cfc ; padding: 10px; border: 1px solid green;">
495 The write-down as one big matrix as above works if and only if the
496 chosen Riemann solver yield a linear combination of the @f$ u^-/u^+ @f$
497 or @f$ x^-/x^+ @f$ values respectively. In many sophisticated solvers,
498 it will not be linear. In this case, the system becomes a non-linear
504 The data storage within a cell is not prescribed by Peano, i.e. you can,
505 in principle decide if you want to reorder the @ref page_multigrid_terminology "nodes within the cell"
506 or you can decide to order the degrees of freedom overall as SoA as opposed
507 to AoS. At the end of the day, Peano maintains an array of doubles per
508 cell and does not really care about the semantics of those. However, we
509 strongly recommend that you stick to the default configuration where
512 @image html degree_of_freedom_semantics.png
514 - Nodes are enumerated following our @ref page_multigrid_terminology "dimension-generic conventions".
515 - Data are stored as an Array of Structs (AoS), i.e. first all the unknowns
516 of the first node, then all of the second, and so forth.
518 If you change any ordering, you have to take this into account for all
519 associated operators (such as projections onto the faces), but you also
520 will have to amend the plotters which are tied to these conventions.
524 - stored as AoS by default;
525 - per face, we first store the left projection @f$ x^- @f$, then the right projection @f$ x^+ @f$
526 and then the outcome of the flux or associated reconstruction. That is, all the @f$ x^- @f$
527 are held en bloc first before we hold @f$ x^+ @f$.
528 - Per @f$ x^\pm @f$, the storage of the nodes and unknowns
529 follows again the @ref page_multigrid_terminology "dimension-generic conventions".
531 It is important to recognise that we store three types of quantities
532 per unknown: Its left projection, its right projection, and its
533 outcome which can be either the flux or a reconstructed
534 income. The semantics of the outcome depend upon your operators, i.e.
535 they are not baked into the formalism. Furthermore, there is no
536 assumption about the spacing along the faces. Naturally, you might
537 assume that a Gauss-Lobatto layout within the cell implies a
538 Gauss-Lobatto distribution on the faces. However, no such thing
539 is encoded within the code, as we solely work with operators. They
540 have to encode the correct spacing.
542 Another assumptions is baked into the realisation as a default:
543 The polynomial degree on the facet result @f$ x_f @f$ equals the
544 polynomial degree within the cell.
546 @todo Different polynomial degrees are simple to realise. This is
547 something to do in the future. Note that +/- are to be
548 taken from the same space as the cell, as they are projections,
549 but maybe we want to alter this, and even provide completely
550 different degrees for the result of the Riemann solve.
556 Multiple examples are shipped with Peano which demonstrate how to use the
557 plain Discontiuous Galerking solver:
559 - benchmarks/multigrid/petsc/poisson/dg.py Hosts a simple solver for the
560 Poisson equation. All the rationale are discussed in @ref benchmarks_documentation_multigrid_petsc_poisson.
573 cell_cell_lhs_matrix,
574 cell_cell_rhs_matrix,
575 cell_cell_lhs_matrix_scaling,
576 cell_cell_rhs_matrix_scaling,
579 cell_to_face_matrix_scaling,
580 face_to_cell_matrix_scaling,
581 face_face_Riemann_problem_matrix,
582 quadrature_points_in_unit_interval,
583 gaussian_integration_points,
584 gaussian_integration_weights
588 Constructor of DG class
591## Cell to face projection
593Consult @ref page_multigrid_terminology "Peano's enumeration for multigrid" page
594for info on the ordering. For 2d with a layout as below
596@image documentation/Doxygen/Multigrid/face-ordering.png
598the matrix is applied to
603 x_0^+ \\ x_1^+ \\ x_2^- \\ x_3^-
613@param cell_cell_lhs_matrix: [Float] or []
614 Pass in [] if you prefer to assemble the matrix per cell yourself. This
615 is the @f$ \tilde{\mathcal{L}} @f$ part above. The ordering of the
616 degrees of freedom follows the discussion above. The matrix is from
617 @f$ K(p+1) \times K(p+1) @f$.
619@param cell_cell_rhs_matrix: [Float] or []
620 Pass in [] if you prefer to assemble the matrix per cell yourself. This
621 matrix arises from @f$ \int f \cdot v \ dx @f$ and is applied to the nodal
622 values of the PDE within the quadrature nodes of a cell. The matrix is
623 from @f$ K(p+1) \times K(p+1) @f$.
625@param cell_cell_lhs_matrix_scaling: Positive integer
626 The lhs matrix is scaled with @f$ h^x @f$ where the x is defined by this
627 parameter. If you don't specify an lhs matrix, i.e. you prefer to inject
628 this matrix manually, you can leave this parameter None, as it is not
631@param cell_cell_rhs_matrix_scaling: Positive integer
632 The rhs matrix is scaled with @f$ h^x @f$ where the x is defined by this
633 parameter. If you don't specify a rhs matrix, i.e. you prefer to inject
634 this matrix manually, you can leave this parameter None, as it is not
637@param cell_to_face_matrix: [Float]
638 If you study the global equation system above, you recognise that
639 this matrix is applied to the left or right solutions of a cell's
642@param face_face_Riemann_problem_matrix: [Float] or []
643 Pass in [] if you prefer to assemble the matrix per cell yourself. This
644 matrix represents the actual Riemann solver. It accepts the left projection
645 @f$ u^- @f$ and @f$ u^+ @f$ and yields a "reconstructed" value at this
646 point passed into the differential operator. Therefore, this is a band
651 super( DiscontinuousGalerkinDiscretisation, self ).
__init__( name,
656 self.polynomial_degree = polynomial_degree
657 self.dimensions = dimensions
658 self.cell_unknowns = cell_unknowns
659 self.face_unknowns = face_unknowns
661 self.cell_cell_lhs_matrix = cell_cell_lhs_matrix
662 self.cell_cell_rhs_matrix = cell_cell_rhs_matrix
663 self.cell_to_face_matrix = cell_to_face_matrix
664 self.face_to_cell_matrix = face_to_cell_matrix
665 self.face_face_Riemann_problem_matrix = face_face_Riemann_problem_matrix
667 self.cell_cell_lhs_matrix_scaling = cell_cell_lhs_matrix_scaling
668 self.cell_cell_rhs_matrix_scaling = cell_cell_rhs_matrix_scaling
669 self.cell_to_face_matrix_scaling = cell_to_face_matrix_scaling
670 self.face_to_cell_matrix_scaling = face_to_cell_matrix_scaling
672 self.quadrature_points_in_unit_interval = quadrature_points_in_unit_interval
673 self.gaussian_integration_points = gaussian_integration_points
674 self.gaussian_integration_weights = gaussian_integration_weights
676 assert( len( quadrature_points_in_unit_interval ) == len( gaussian_integration_points ) == len( gaussian_integration_weights ),
"Something has gone wrong. These should all be same length for fixed polynomial degree" )
687 def nodes_per_cell(self):
688 return (self.polynomial_degree + 1) ** self.dimensions
692 def nodes_per_face(self):
693 return (self.polynomial_degree + 1) ** (self.dimensions-1)
696 def add_to_Peano4_datamodel( self, datamodel, verbose ):
699 Initialise index model
701 There is an index data model, where all the vertices, faces and cells
702 are properly enumerated. This one is taken care of by the superclass.
703 The thing we have to add our our PDE-specific unknowns.
706 super( DiscontinuousGalerkinDiscretisation, self ).add_to_Peano4_datamodel(datamodel,verbose)
707 datamodel.add_cell(self._cell_pde_data)
710 def add_use_statements(self, observer):
711 super( DiscontinuousGalerkinDiscretisation, self ).add_use_statements(observer)
712 observer.use_cell(self._cell_pde_data)
716 def add_to_plot(self, observer):
719 Tell the project's observer how to plot the data
721 Nothing fancy here. We add plotters from Peano's toolbox to visualise
722 solution and right-hand side.
725 observer.add_action_set( PlotDGDataInPeanoBlockFormat(self) )
729 def add_to_create_grid(self, observer):
734 def add_to_enumerate_and_init_solution(self, observer):
737 Solution initialisation
739 Close to trivial: Just add the action set petsc.actionsets.InitVertexDoFs
743 observer.add_action_set( EnumerateDoFs(self,
744 boundary_vertices_hold_data =
False,
745 boundary_faces_hold_data =
True,
749 observer.add_action_set( InitCellDoFs(self) )
752 def add_to_assemble(self, observer):
754 observer.add_action_set( ImposeDirichletBoundaryConditionsWithInteriorPenaltyMethod(self) )
757 def add_to_init_petsc(self, observer):
760 def add_to_map_solution_onto_mesh(self, observer):
761 observer.add_action_set( ProjectPETScSolutionOnCellsBackOntoMesh(self) )
762 observer.add_action_set( PlotExactSolution(self) )
765 def add_implementation_files_to_project(self, namespace, output, subdirectory=""):
768 The solver creates two classes: An abstract base class and its
769 implementations. The abstract base class will be overwritten if
770 there is one in the directory. I pipe all the Python constants
771 into this base class, so they are available in the C++ code.
773 The implementation file will not be overwritten, as I assume that
774 the users will write their own domain-specific stuff in there. If
775 it does not exist yet, I'll create an empty stub which can be
776 befilled with something meaningful.
778 Besides the creation of these two files, I also add the files to the
779 project artifacts and the makefile, so they are captured by the build
783 templatefile_prefix = os.path.dirname( os.path.realpath(__file__) ) +
"/"
784 templatefile_prefix += self.__class__.__name__
790 "SOLVER_INCLUDES":
"",
791 "MIN_H" : self.min_h,
792 "MAX_H" : self.max_h,
793 "POLYNOMIAL_DEGREE": self.polynomial_degree,
794 "CELL_UNKNOWNS": self.cell_unknowns,
795 "FACE_UNKNOWNS": self.face_unknowns,
796 "CELL_CELL_LHS_MATRIX": self.cell_cell_lhs_matrix,
797 "CELL_CELL_RHS_MATRIX": self.cell_cell_rhs_matrix,
798 "FACE_FACE_RIEMANN_PROBLEM_MATRIX": self.face_face_Riemann_problem_matrix,
799 "CELL_TO_FACE_MATRIX": self.cell_to_face_matrix,
800 "FACE_TO_CELL_MATRIX": self.face_to_cell_matrix,
801 "QUADRATURE_POINTS_IN_UNIT_INTERVAL": self.quadrature_points_in_unit_interval,
802 "GAUSSIAN_INTEGRATION_POINTS" : self.gaussian_integration_points,
803 "GAUSSIAN_INTEGRATION_WEIGHTS" : self.gaussian_integration_weights,
804 "CELL_CELL_LHS_MATRIX_SCALING": self.cell_cell_lhs_matrix_scaling,
805 "CELL_CELL_RHS_MATRIX_SCALING": self.cell_cell_rhs_matrix_scaling,
806 "CELL_TO_FACE_MATRIX_SCALING": self.cell_to_face_matrix_scaling,
807 "FACE_TO_CELL_MATRIX_SCALING": self.face_to_cell_matrix_scaling,
808 "SOLVER_NAME": self.typename(),
812 templatefile_prefix +
".Abstract.template.h",
813 templatefile_prefix +
".Abstract.template.cpp",
814 "Abstract" + self.typename(),
821 templatefile_prefix +
".template.h",
822 templatefile_prefix +
".template.cpp",
830 output.add( generated_abstract_header_file )
831 output.add( generated_solver_files )
832 output.makefile.add_h_file( subdirectory +
"Abstract" + self._name +
".h", generated=
True )
833 output.makefile.add_h_file( subdirectory + self._name +
".h", generated=
True )
834 output.makefile.add_cpp_file( subdirectory +
"Abstract" + self._name +
".cpp", generated=
True )
835 output.makefile.add_cpp_file( subdirectory + self._name +
".cpp", generated=
True )
839 def number_of_matrix_entries_per_vertex(self):
842 Nothing is associated with the vertex. There's nothing to be done here.
848 def number_of_matrix_entries_per_face(self):
851 In the DG scheme, we have a projection from the left and we have a
852 projection from the right. These values are proper projections, i.e.
853 they do not carry any semantics. Then we have to solve the Riemann
854 problem, and need one more unknown to store the outcome of that one.
857 return 3 * self.face_unknowns * self.nodes_per_face
860 def number_of_matrix_entries_per_cell(self):
863 This is the actual unknowns per cell
866 return self.cell_unknowns * self.nodes_per_cell
Trigger assembly of PETSc matrix and also assemble rhs vector.
get_constructor_body(self)
Define body of constructor.
get_action_set_name(self)
Configure name of generated C++ action set.
get_attributes(self)
Return attributes as copied and pasted into the generated class.
str TemplateTouchCellFirstTime
str TemplateTouchFaceFirstTime
get_includes(self)
Consult petsc.Project for details.
__init__(self, solver)
Yet to be written.
get_body_of_operation(self, operation_name)
Provide C++ code snippet for peano4.solversteps.ActionSet.OPERATION_TOUCH_VERTEX_FIRST_TIME Provide...
user_should_modify_template(self)
The action set that Peano will generate that corresponds to this class should not be modified by user...
Abstract base class for all PETSc solvers.
Specialisation of dastgen2.attributes.DoubleArray which relies on Peano's tarch.
Default superclass for any data model in Peano which is stored within the grid.
Action set (reactions to events)