 |
Peano
|
|
Peano
|
This tutorial runs through the basic steps how to construct a geometric single level Jacobi solver for
\( -\Delta u = f \)
While noone would use such a numerical scheme as a real solver for any problem of relevant size, this type of solver often is a building block of proper multigrid algorithms.
MGHyPE codes need to know where to find the MGHyPE and Peano Python sources. For this, we can either add the Python sources to our environment, or export the path variable prior to invoking our script/kicking off our Jupyter notebook, or we set the path hard-coded in our script. To be honest, this last approach is likely the least elegant one, as it hard-codes the pathes, but it works:
Once we have ensured that all pathes are correct, we import the packages of interest:
Every MGHyPE project starts with an MGHyPE project object:
This creates the domain-specific project. We immediately set the domain properties:
We parse the system setup next, so all the compiler settings and stuff is known to MGHyPE/Peano. As we develop code, we switch to asserts mode, which means a lot of checks are injected into the code base. Once we are sure everything works, we toggle this switch to Release:
Once we have a project, we can start to construct our solver. One of the first things that we will need is the element assembly matrices. MGHyPE always assumes that the discretisation yields a system
\( Ax = Mb \)
with a system matrix \( A\) and a right-hand side matrix \( M\) . If we employ the scheme as an actual solver, this \( M\) is the mass matrix corresponding to an integral
\( \int (f, \varphi ) \, dx \)
Furthermore, both \( A\) and \( M\) are always handed over to the solver as normalised matrices (over the unit cube) plus some h-scaling. There's no need to create these matrices manually. We can use one of MGHyPE's factory mechanisms:
We briefly inspect the results (something we would likely not do for a "real" simulation):
The first matrix is the standard (2d) element matrix for the Poisson equation. The scaling array indicates that this one has to be scaled with \( h^0 \) before we continue. More details can be found in the documentation of these routines on the webpage. When we study the mass matrix, there are no big surprises either.
Next we create the solver using these matrices and add it to the project:
We have now set up the solver. Hence, we can lower it to C++:
We build the project now:
So far, this is a dysfunctional solver, as we haven't set any boundary conditions, and no right-hand side. There are multiple ways to do this in MGHyPE. We demonstrate the low-level version here: Every solver is mapped onto an AbstractCollocatedPoisson class and a subclass called CollocatedPoisson. The abstract class is regenerated every time we run Python and basically stores all the data from the Python within the C++ code. It represents a translation from Python into the compiled world. The subclass can be used by us to tailor the code to our actual needs.
Open the C++ file of the solver and ensure that it realises a proper solver, i.e. something alike
We next run the experiment but ensure that the output is piped into a file, as we'll later on use this dump to run some analysis:
Nothing controversial here. The residual goes down albeit slowly, and we therefore never reach our numerical termination criterion. Instead, we quit once we exceed the maximum number of iterations.
Peano writes its own bespoke data file format. There are some scripts delivered to convert them into VTK which you can load into Paraview or VisIt, e.g. we use exactly this script here to convert the outcome into a standard visualisation format. Unfortunately, my pvpython sometimes seems to miss out on environment variables. So, I have to re-add those of interest manually:
We can now visualise the pvd/vtu files in Paraview or VisIt. As this is all only 2d, I use the PyVista package to do this directly within the Jupyter notebook. The first step grabs the package and therefore is required once and once only. The second one then creates the pics:
No surprises with the right-hand side. We however see that the result is nowhere near the real solution. This is due to the poor convergence rates. In the script below, we might want to pass in clim=[from,to] to limit the range displayed. Note furthermore that we pick a second time step here, and do not load the pvd file. Each solver, by default, writes the initial solution and the last one - unless we make it plot every single mesh sweep.
We close this first tutorial with a plot of the convergence behaviour. For this, we open the output that we'd streamed into a text file and parse the results. As the dumps are kind of standardised, we can use some utility functions shipped with MGHyPE to do the actual parsing:
1.9.8