dg.py

Go to the documentation of this file.

import os, sys
import argparse
import peano4
import mghype
import numpy as np
 
 
def get_args():
  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("-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 (default 2)")
  parser.add_argument("-p",   "--poly_degree",          dest="poly_degree",    type=int, default=1, help="Polynomial Degree (default 1)")
  parser.add_argument("-meshsize",  "--meshsize",       dest="meshsize",       default=0.3, help="Mesh size (default 0.3)")
  parser.add_argument("-omega_c",  "--omega_c",         dest="omega_c",        default=0.5, help="Mesh size (default 0.5)")
  parser.add_argument("-omega_f",  "--omega_f",         dest="omega_f",        default=1.0, help="Mesh size (default 1.0)")
  parser.add_argument("-m",   "--mode",                 dest="mode",           choices=["release","stats","asserts", "trace"], required=True, help="Pick build type" )
  parser.add_argument("-t",   "--test-run",             dest="test_run",       action="store_true", default=False, help="Trigger test run")
  
  return parser.parse_args()
 
      
def configure_and_build_project(args,project):
  #
  # Configure build
  #
  if args.mode=="release":
    build_mode = peano4.output.CompileMode.Release
  if args.mode=="stats":
    build_mode = peano4.output.CompileMode.Stats
  if args.mode=="asserts":
    build_mode = peano4.output.CompileMode.Asserts
  if args.mode=="trace":
    build_mode = peano4.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                   = 10000,
#    max_iterations                   = 10,
  )
  
  project.set_load_balancing( "toolbox::loadbalancing::strategies::SpreadOutOnceGridStagnates", "new ::exahype2::LoadBalancingConfiguration(0.9, 0, false, 8, 50)" )
  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()
 
  
def run_and_validate(args,executable):
  os.system( "make clean" )
  if args.j>0:
    os.system( "make -j{}".format(args.j) )
  else:
    os.system( "make -j$(nproc --all)" )
    
  if args.test_run:
    os.system( "./" + executable + " > " + executable + ".out")
    output_file = open(executable+".out", "r")
    initial_r_max = 0.0
    final_r_max   = 0.0
    for line in output_file.readlines():
      if "endMeshSweep()" in line and initial_r_max==0.0:
        initial_r_max = float( line.split( "|r|_max=" )[1].split(",")[0] )
      if "endMeshSweep()" in line:
        final_r_max = float( line.split( "|r|_max=" )[1].split(",")[0] )
    print( "initial residual = {}".format( initial_r_max) )
    print( "final residual = {}".format( final_r_max) )
    if initial_r_max>0.0:
        improvement = initial_r_max / final_r_max;
    else: 
        improvement = 0.0
    print( "improvement = {}".format( improvement) )
    if improvement>10.0:
      print( "SUCCESS")
    else:
      sys.exit( "FAILED: residual not within bounds" )
 
 
 
if __name__ == "__main__":
    args = get_args()
    
    
    executable = "peano4matrixfree"
        
    project = mghype.matrixfree.api.Project( project_name = "DGPointwise", 
                                  namespace = [ "tests", "multigrid", "matrixfree", "poisson" ],
                                  executable = executable
                                )
    
    
    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()
    
    
    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 = 0.01,
      integration_points = 0.5*(matrices.get_nodal_points_1d() + 1.0) # Should be specified in [0,1] interval
    )
    
    
    project.add_solver(solver)
    
    configure_and_build_project(args,project)
    run_and_validate(args,executable)