 |
Peano
|
|
Peano
|
This tutorial runs through the steps that we perform to write a single-level DG solver on a regular grid. We implement the standard "benchmark" for this type of solvers, which is solving the Poisson equation on the unit square. All sources for this can be found in tests/multigrid/matrixfree/poisson/dg.py, as the solver is routinely used to do convergence measurements. While we follow some steps of this script, the test Python script has more technical flavours/details which we omit in the notebook. Therefore, it is reasonable to study the present text in combination with the actual script once you've done the notebook once.
The notebook clarifies how to construct an MGHyPE project with one solver of the type multigrid.matrixfree.solvers.api.DiscontinuousGalerkinDiscretisationPointWiseRiemannSolver (cmp multigrid/matrixfree/solvers/api/DiscontinuousGalerkinDiscretisationPointWiseRiemannSolver.py in src).
Since we solve the Poisson equation, and employ a Discontinuous Galerkin formalism with a point-wise Riemann solver, we start from
\( \mathcal{L}(u) = - \Delta u = f \)
and apply a weak formulation:
\( \begin{eqnarray} \int \left( - \Delta u, v \right) dx & = & \int (f, v) dx \\ \int _c \left( \nabla u, \nabla v \right) dx - \sum _f \int _f \left( \nabla u, n \right) v dS(x) & = & \int (f, v) dx \\ \int _c \left( \nabla u, \nabla v \right) dx - \sum _f \int _f \left( \mathcal{R}(u^-,u^+), n \right) \phi dS(x) & = & \int (f,v) dx. \end{eqnarray} \)
Before we start, we assume that you have either set export PYTHONPATH=../../../../python:../../../../src or you invoke
We create an MGHyPE project (bear in mind that we basically discuss how to write the test case in the repo; we'd likely use a different namespace otherwise):
Peano 4 (C) www.peano-framework.org VTK is not available, not loading peano4.visualisation. ParaView is not available, not loading peano4.visualisation. Peano 4 - Multigrid extension project.subdirectories = ['']
The solver construction also is a "one-liner" using the following observations:
Warning: penalty parameter gamma_F is set to 2/h and does not depend on the polynomial degree!
This code block employs MGHyPE's matrix generator quite heavily. It makes sense to consult this class' documentation for furhter info. Or switch from Gauss Lobatto (GL) to Gausss Legendre.
@Alex can you add a remark how to do so?
We next add this solver to the project and set some global parameters:
In line with the hitchhiker's guide, we now lower this into a Peano project and make the latter dump the C++ code:
parse configure outcome ../../../config.log to extract configure settings
parse configure outcome ../../../src/Makefile to extract compile settings
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[6], line 1
----> 1 peano4_project = project.generate_Peano4_project(False)
2 peano4_project.output.makefile.add_cpp_file( "Scenario.cpp" )
3 peano4_project.generate()
File ~/git/Peano/src/mghype/matrixfree/api/Project.py:176, in Project.generate_Peano4_project(self, verbose)
173 self._project.solversteps.add_step(self.algorithm_step_plot)
175 self.__configure_makefile()
--> 176 self.__generate_solver_repository()
178 self._project.output.readme.add_package_description( """
179
180
181 """ )
183 return self._project
File ~/git/Peano/src/mghype/matrixfree/api/Project.py:204, in Project.__generate_solver_repository(self)
196 templatefile_prefix = os.path.realpath(__file__).replace( ".pyc", "" ).replace( ".py", "" )
197 template_name = templatefile_prefix + ".SolverRepository.template"
198 generated_repository_files = peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
199 template_name + ".h",
200 template_name + ".cpp",
201 "SolverRepository",
202 self._project.namespace, # + ["repositories"],
203 self._project.subdirectory + "repositories",
--> 204 self.dictionary,
205 True,
206 True)
208 self._project.output.add( generated_repository_files )
209 self._project.output.makefile.add_h_file( self._project.subdirectory + "repositories/SolverRepository.h", generated=True )
AttributeError: 'Project' object has no attribute 'dictionary'
The only "new" thing here's the line add_cpp_file() which we discuss next.
Whenever we create a solver, MGHyPE creates two classes in our home directory. As we strictly split each class into its header and its implementation, we get four files:
01_matrix-free-cg.dox 02_matrix-free-dg.ipynb solver-class-hierarchy.png 01_matrix-free-cg.ipynb 11_petsc-cg.dox solver-class-hierarchy.svg 02_matrix-free-dg.dox 11_petsc-cg.ipynb
The abstract class inherits from some class with MGHyPE. It will contain all the matrices and whatever stuff we have specified in the Python script. The class Poisson.cpp is the one where we can inject our own behaviour.
For example, we should use it to implement our own initial guess and the right-hand side. See
for an example. We can always add additional files to the Peano project, and we've done exactly this through the add_cpp_file() above. This allows us to use the same right-hand side, i.e. the same scenario, in each and every test case.
Open the generated .h and .cpp file of the implementation, compare its content to the files in the test directory, and implement your own solver.
The compilation is straightforward. We simply invoke the Makefile. Once the implementation of the solver is in place, we can rerun the Python script arbitrarily often. The actual implementation will never ever been overwritten, i.e. we can play around with different polynomial degrees, matrices, and so forth, within the Python script.
If the make is successful, we will obtain an executable, which we can once again call from the Jupyter notebook. Note that I pipe the output into a file, so I can easily access it later on:
01_matrix-free-cg.dox 02_matrix-free-dg.ipynb solver-class-hierarchy.png 01_matrix-free-cg.ipynb 11_petsc-cg.dox solver-class-hierarchy.svg 02_matrix-free-dg.dox 11_petsc-cg.ipynb /bin/bash: line 1: ./peano: No such file or directory
MGHyPE comes along with some scripts to parse the dump of a solver and to visualise it directly via matplotlib. As we have piped the solver's output into a file output.txt, we can directly invoke the script from within the Jupyter notebook:
Obviously, visualising the convergence behaviour is only part of the game. We also want to see what the solution looks like. For this, we have to convert the output of the program. Peano provides some scripts for this. It is tricky to embed them into a Jupyter notebook, so they are best invoked on the command line:
From hereon, you can open Paraview or VisIt and directly visualise the outcome. You should obtain a picture similar to
1.10.0