This documentation is autogenerated from ./tutorials/mghype/examples/02_higher-order-cg.ipynb. Once you clone the Peano repository, please work against this notebook or study Peano's generic description how to generate web docu from Jupyter notebooks.

We stick to

\( -\Delta u = f \)

but solve it with a higher order discretisation which is subject to a block Jacobi: All the unknowns within a cell are updated via direct inversion, while the unknowns on the facets are still updated by a Jacobi scheme.

Prepare the environment

We kick off wiht some system settings:

import sys, os
 
sys.path.insert(0, os.path.abspath("../../../python/"))
sys.path.insert(0, os.path.abspath("../../../src/"))

and then start the actual project:

import peano4.api, mghype
 
project = mghype.matrixfree.api.Project( project_name = "Tutorial02", 
                                         namespace = [ "tutorials", "mghype", "examples" ],
                                         executable = "tutorial02",
                                        )
cube_size = 1.0
dimensions = 2
project.set_global_simulation_parameters(
  dimensions            = dimensions,
  offset                = [ 0.0       for _ in range(dimensions) ],
  domain_size           = [ cube_size for _ in range(dimensions) ],
  plot_after_each_mesh_sweep = True,
  max_iterations             = 20
)
 
project.set_Peano4_installation( "../../..", peano4.api.output.CompileMode.Asserts )

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:

Add a solver

We follow the discussion on Continuous Galerkin for the Poisson equation and create a solver for a scalar system, which hosts the residual and the diagonal on the faces, i.e. two unknowns there:

polynomial_degree = 3
 
poisson_equation_matrices = mghype.api.matrixgenerators.GLMatrixFree(
  dimensions  = dimensions,
  poly_degree = polynomial_degree,
  unknowns_per_cell_node = 1,
  solutions_per_face_node = 1,   # only the weight of the shape function
  projections_per_face_node = 1  # residual plus diagonal value
)
 
system_matrix, system_matrix_scaling = poisson_equation_matrices.get_cell_system_matrix_for_laplacian()
mass_matrix,   mass_matrix_scaling   = poisson_equation_matrices.get_cell_mass_matrix()
#projection_matrix, projection_matrix_scaling = poisson_equation_matrices.get_face_from_cell_matrix()

Again, we can next create our solver:

solver_name = "CGHigherOrderPoisson"
mesh_size = 0.1
solver = mghype.matrixfree.solvers.api.ContinuousGalerkinDiscretisationGaussLobattoWithBlockJacobi(solver_name,
                                                              1, # unknowns per node, i.e. scalar equation
                                                              dimensions,
                                                              polynomial_degree,
                                                              mesh_size,
                                                              mesh_size,
                                                              assembly_matrix = system_matrix,
                                                              assembly_matrix_scaling=system_matrix_scaling,
                                                              mass_matrix = mass_matrix,
                                                              mass_matrix_scaling = mass_matrix_scaling,
                                                              cell_relaxation = 0.6,                      
                                                              solver_tolerance=0.01,
                                                              integration_points=poisson_equation_matrices.get_integration_points_over_unit_interval()
                                                            )

project.add_solver(solver)

Create and build project

We have now set up the solver. Hence, we can lower it to C++:

peano4_project = project.generate_Peano4_project(verbose=True)

peano4_project.generate()

We build the project now:

!make -j8

So far, this is a disfunctional 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.

Exercise

Open the C++ file of the solver and ensure that it realises a proper solver, i.e. something alike

#include <cmath>
 
 
tutorials::mghype::examples::CollocatedPoisson::CollocatedPoisson() {
}
 
 
tutorials::mghype::examples::CollocatedPoisson::~CollocatedPoisson() {
}
 
 
void tutorials::mghype::examples::CollocatedPoisson::initVertex(
  const tarch::la::Vector<Dimensions, double>&  x,
  const tarch::la::Vector<Dimensions, double>&  h,
  tarch::la::Vector< 1, double >&  value,
  tarch::la::Vector< 1, double >&  rhs
) {
  value = 0.0;
  rhs   = 1.0;
  for (int d=0; d<Dimensions; d++) {
    rhs *= std::sin( tarch::la::PI * x(d) );
  } 
}
 
void tutorials::mghype::examples::CollocatedPoisson::setBoundaryConditions(
  const tarch::la::Vector<Dimensions, double>&  x,
  const tarch::la::Vector<Dimensions, double>&  h,
  tarch::la::Vector< 1, double >&  value
) {
  value = 0.0;
}

Running the experiment

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:

!rm *.peano-patch-file
!./tutorial02 > tutorial02.out
!cat tutorial02.out

Nothing controversial here. The residual goes down albeit slowly, and we therefore never reach out numerical termination criterion. Instead, we quit once we exceed the maximum number of iterations.

Visualisation

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:

import subprocess
import os
my_env = os.environ
my_env["PYTHONPATH"] = "../../../python:" + my_env["PYTHONPATH"]
subprocess.call(["pvpython", "../../../python/peano4/visualisation/render.py", "CGHigherOrderPoisson.rhs.peano-patch-file"], env=my_env)
subprocess.call(["pvpython", "../../../python/peano4/visualisation/render.py", "CGHigherOrderPoisson.solution.peano-patch-file"], env=my_env)

Display the outcome

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:

import sys

!{sys.executable} -m pip install --break-system-packages pyvista[all]

import pyvista
mesh = pyvista.read( "CGHigherOrderPoisson.rhs.pvd" )
mesh.plot()

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.

plotter = pyvista.Plotter()
mesh = pyvista.read( "CGHigherOrderPoisson.solution.peano-patch-file-1.vtu" )
plotter.add_mesh(mesh, show_edges=True ) # , clim=[0.0, 0.01]
plotter.show()

Study the convergence behaviour

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:

res_max, res_2, res_h = mghype.api.postprocessing.parse_residuals( "tutorial02.out", solver_name )

du_max, du_2, du_h = mghype.api.postprocessing.parse_solution_updates( "tutorial02.out", solver_name )

import matplotlib.pyplot as plt
plt.plot(res_2,   label='$|res|_2$',       marker='<')
plt.plot(du_2,    label='$|\delta u|_2$',  marker='>')
#plt.plot(res_max, label='$|res|_{max}$', marker='>')
#plt.plot(res_h,   label='$|res|_h$',     marker='v')
plt.yscale('log')
plt.xlabel('iterations')
plt.ylabel('residual')
plt.legend()
plt.savefig( "tutorial02.pdf" )
plt.savefig( "tutorial02.png" )
plt.show()