dg-two-grid-h-coarsen.py

Go to the documentation of this file.

import os, sys
import argparse
import peano4.api
import mghype
import numpy as np
 
print( """
 
Peano 4 multigrid solver using Discontinuous Galerkin
 
""")
 
parser = argparse.ArgumentParser(description='Multigrid - Poisson solver')
parser.add_argument("-j",   "--parallel-builds",    dest="j",                  type=int, default=-1, help="Parallel builds" )
parser.add_argument("-pt",   "--plot-sweep",     dest="plot_after_each_mesh_sweep", action="store_true", default=False, help="Plot after each timestep. By default (False) we plot after initialisation and at the end" )
parser.add_argument("-fas",   "--use-fas",          dest="fas", action="store_true", default=False, help="Switch to FAS" )
parser.add_argument("-pd",  "--peano-dir",          dest="peanodir",       default="../../../../", help="Peano4 directory" )
parser.add_argument("-v",   "--verbose",            dest="verbose",        action="store_true", default=False, help="Verbose")
parser.add_argument("-d",   "--dimensions",         dest="dimensions",     type=int, default=2, help="Dimensions")
parser.add_argument("-p",   "--poly_degree",        dest="poly_degree",    type=int, default=2, help="Polynomial Degree")
parser.add_argument("-meshsize",  "--meshsize",     dest="meshsize",       default=0.3, help="Mesh size")
parser.add_argument("-omega_c",  "--omega_c",       dest="omega_c",       default=0.5, help="Smoothing parameter for cells")
parser.add_argument("-omega_f",  "--omega_f",       dest="omega_f",       default=1.0, help="Smoothing parameter for faces")
parser.add_argument("-m",   "--mode",               dest="mode",           choices=["release","stats","asserts", "trace"], required=True, help="Pick build type" )
args = parser.parse_args()
 
project = mghype.matrixfree.api.Project( project_name = "DGTwoGrid", 
                              namespace = [ "tests", "multigrid", "matrixfree", "poisson" ]
                            )
 
matrices = mghype.api.matrixgenerators.GaussLobattoMatrixFree(
  args.dimensions,
  args.poly_degree,
  1,  # Unknowns per cell dof. Scalar PDE here
  2,  # We use the penalty formulation (see docu in tutorials)
  2
)
 
assembly_matrix, assembly_matrix_scaling = matrices.get_cell_system_matrix_for_laplacian()
mass_matrix, mass_matrix_scaling         = matrices.get_cell_mass_matrix()
face_from_cell_projection, \
face_from_cell_projection_scaling        = matrices.get_face_from_cell_matrix()
cell_from_face_projection, \
cell_from_face_projection_scaling        = matrices.get_cell_from_face_matrix()
approximate_system_matrix, \
approximate_system_matrix_scaling        = matrices.get_A_tilde()
 
#! [Instantiate solver]
solver = mghype.matrixfree.solvers.api.DiscontinuousGalerkinDiscretisationPointWiseRiemannSolver(
  name = "DGPoisson",
  dimensions = args.dimensions,
  poly_degree = args.poly_degree,
  unknowns_per_cell_node = 1,
  solutions_per_face_node = 2,
  projections_per_face_node = 2,
  min_h = args.meshsize,
  max_h = args.meshsize,
  assembly_matrix= assembly_matrix, 
  assembly_matrix_scaling = assembly_matrix_scaling, 
  mass_matrix = mass_matrix, 
  mass_matrix_scaling = mass_matrix_scaling, 
  face_from_cell_projection = face_from_cell_projection, 
  face_from_cell_projection_scaling = face_from_cell_projection_scaling,
  cell_from_face_projection = cell_from_face_projection, 
  cell_from_face_projection_scaling = cell_from_face_projection_scaling,
  riemann_matrix = matrices.get_face_face_riemann_problem_matrix(),
  boundary_matrix = matrices.get_boundary_matrix(),
  cell_relaxation = args.omega_c,
  face_relaxation = args.omega_f,  # default here is 1, i.e. we solve exactly. See class docu
  approximate_system_matrix = approximate_system_matrix,
  approximate_system_matrix_scaling = approximate_system_matrix_scaling,
  solver_tolerance = 1.0e-6,
  integration_points = 0.5*(matrices.get_nodal_points_1d() + 1.0) # Should be specified in [0,1] interval
)
#! [Instantiate solver]
 
#! [Construct correction matrices]
correction_equation_matrices = mghype.api.matrixgenerators.DLinearMassIdentity(
  dimensions = args.dimensions,
  poly_degree = 1,
  unknowns_per_node = 1
)
 
intergridOperators = mghype.api.matrixgenerators.blockmatrix.IntergridOperatorsH()
 
system_matrix, system_matrix_scaling = correction_equation_matrices.get_cell_system_matrix_for_laplacian()
mass_matrix,   mass_matrix_scaling   = correction_equation_matrices.get_cell_mass_matrix()
 
correction_matrix, correction_matrix_scaling = correction_equation_matrices.get_cell_mass_matrix()
#! [Construct correction matrices]
 
"""
We use the identity matrix everywhere...
"""
correction_equation_solver = mghype.matrixfree.solvers.api.CollocatedMGSingleThreaded("CollocatedMG",
                                                    1,
                                                    args.dimensions,
                                                    args.meshsize,
                                                    args.meshsize,
                                                    local_assembly_matrix = system_matrix,
                                                    local_assembly_matrix_scaling=system_matrix_scaling,
                                                    mass_matrix = mass_matrix,
                                                    mass_matrix_scaling = mass_matrix_scaling,
                                                    correction_matrix = correction_matrix,
                                                    correction_matrix_scaling = correction_matrix_scaling,
                                                    solver_tolerance=1e-12,
                                                    pre_smoothing_steps =6,
                                                    post_smoothing_steps=5,
                                                    restriction_matrix=intergridOperators.restriction(),
                                                    prolongation_matrix=intergridOperators.prolongation(),
                                                    descend_order=0,
                                                    )
 
#! [Coupling]
intergrid_operators = mghype.api.matrixgenerators.blockmatrix.IntergridOperators(
  dim = args.dimensions,
  coarse_order = 1,
  fine_order = args.poly_degree
  )
prolongation_matrix = [intergrid_operators.prolongation()]
prolongation_matrix_scaling = [0]
restriction_matrix = [intergrid_operators.restriction()]
restriction_matrix_scaling = [0]
injection_matrix = [intergrid_operators.injection()]
injection_matrix_scaling = [0]
 
solver.preprocessing_action_set = mghype.matrixfree.solvers.api.actionsets.MultiplicativeDGCGCoupling(
    solver,                      # dg_solver
    correction_equation_solver,  # cg_solver
    prolongation_matrix,
    prolongation_matrix_scaling,
    restriction_matrix,
    restriction_matrix_scaling,
    injection_matrix,
    injection_matrix_scaling,
    args.fas
)
#! [Coupling]
 
 
#! [Add solvers]
project.add_solver(correction_equation_solver, auxiliary=True)
project.add_solver(solver)
#! [Add solvers]
 
#
# Configure build
#
if args.mode=="release":
  build_mode = peano4.api.output.CompileMode.RELEASE
if args.mode=="stats":
  build_mode = peano4.api.output.CompileMode.STATS
if args.mode=="asserts":
  build_mode = peano4.api.output.CompileMode.ASSERTS
if args.mode=="trace":
  build_mode = peano4.api.output.CompileMode.TRACE
 
 
cube_size = 1.0
project.set_global_simulation_parameters(
  dimensions            = args.dimensions,
  offset                = [ 0.0       for _ in range(args.dimensions) ],
  domain_size           = [ cube_size for _ in range(args.dimensions) ],
  plot_after_each_mesh_sweep       = args.plot_after_each_mesh_sweep,
  max_iterations = 2000
)
 
project.set_Peano4_installation( args.peanodir, build_mode )
 
peano4_project = project.generate_Peano4_project(args.verbose)
peano4_project.output.makefile.add_cpp_file( "Scenario.cpp" )
peano4_project.generate()
 
os.system( "make clean" )
if args.j>0:
  os.system( "make -j{}".format(args.j) )
else:
  os.system( "make -j" )