started tutorial

Author: pancetta

pySDC/playgrounds/PETSc/playground_pySDC_GS_multiimplicit.py

@@ -1,54 +1,46 @@
 import sys
-
 from mpi4py import MPI
-import time
-from petsc4py import PETSc
-
 import numpy as np
 
-from pySDC.implementations.problem_classes.GrayScott_2D_PETSc_implicit_periodic import petsc_grayscott
-from pySDC.implementations.datatype_classes.petsc_dmda_grid import petsc_data
+from pySDC.implementations.problem_classes.HeatEquation_2D_PETSc_forced import heat2d_petsc_forced
+from pySDC.implementations.datatype_classes.petsc_dmda_grid import petsc_data, rhs_imex_petsc_data
 from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
-from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
 from pySDC.implementations.transfer_classes.TransferPETScDMDA import mesh_to_mesh_petsc_dmda
 from pySDC.implementations.controller_classes.allinclusive_multigrid_MPI import allinclusive_multigrid_MPI
-from pySDC.implementations.controller_classes.allinclusive_multigrid_nonMPI import allinclusive_multigrid_nonMPI
 
 from pySDC.helpers.stats_helper import filter_stats, sort_stats
 
-def main():
 
+def main():
+    """
+    Program to demonstrate usage of PETSc data structures and spatial parallelization,
+    combined with parallelization in time.
+    """
     # set MPI communicator
     comm = MPI.COMM_WORLD
 
     world_rank = comm.Get_rank()
     world_size = comm.Get_size()
 
+    # split world communicator to create space-communicators
     if len(sys.argv) == 2:
         color = int(world_rank / int(sys.argv[1]))
     else:
         color = int(world_rank / 1)
-
     space_comm = comm.Split(color=color)
-    space_rank = space_comm.Get_rank()
-    space_size = space_comm.Get_size()
 
+    # split world communicator to create time-communicators
     if len(sys.argv) == 2:
         color = int(world_rank % int(sys.argv[1]))
     else:
         color = int(world_rank / world_size)
-
     time_comm = comm.Split(color=color)
-    time_rank = time_comm.Get_rank()
-    time_size = time_comm.Get_size()
-
-    print("IDs (world, space, time):  %i / %i -- %i / %i -- %i / %i" % (world_rank, world_size, space_rank, space_size,
-                                                                        time_rank, time_size))
 
     # initialize level parameters
     level_params = dict()
     level_params['restol'] = 1E-08
-    level_params['dt'] = 1.0
+    level_params['dt'] = 0.125
     level_params['nsweeps'] = [1]
 
     # initialize sweeper parameters
@@ -60,51 +52,47 @@ def main():
 
     # initialize problem parameters
     problem_params = dict()
-    problem_params['Du'] = 1.0
-    problem_params['Dv'] = 0.01
-    problem_params['A'] = 0.09
-    problem_params['B'] = 0.086
-    problem_params['nvars'] = [(127, 127)]  # number of degrees of freedom for each level
-    problem_params['comm'] = space_comm
-    problem_params['sol_tol'] = 1E-10
-    problem_params['sol_maxiter'] = 100
+    problem_params['nu'] = 1.0  # diffusion coefficient
+    problem_params['freq'] = 2  # frequency for the test value
+    problem_params['nvars'] = [(129, 129), (65, 65)]  # number of degrees of freedom for each level
+    problem_params['comm'] = space_comm  # pass space-communicator to problem class
+    problem_params['sol_tol'] = 1E-12  # set tolerance to PETSc' linear solver
 
     # initialize step parameters
     step_params = dict()
     step_params['maxiter'] = 50
 
     # initialize space transfer parameters
-    # space_transfer_params = dict()
-    # space_transfer_params['rorder'] = 2
-    # space_transfer_params['iorder'] = 2
-    # space_transfer_params['periodic'] = True
+    space_transfer_params = dict()
+    space_transfer_params['rorder'] = 2
+    space_transfer_params['iorder'] = 2
+    space_transfer_params['periodic'] = False
 
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 20
     # controller_params['predict'] = False
-    # controller_params['hook_class'] = error_output
 
     # fill description dictionary for easy step instantiation
     description = dict()
-    description['problem_class'] = petsc_grayscott # pass problem class
+    description['problem_class'] = heat2d_petsc_forced  # pass problem class
     description['problem_params'] = problem_params  # pass problem parameters
-    description['dtype_u'] = petsc_data  # pass data type for u
-    description['dtype_f'] = petsc_data  # pass data type for f
-    description['sweeper_class'] = generic_implicit  # pass sweeper (see part B)
+    description['dtype_u'] = petsc_data  # pass PETSc data type for u
+    description['dtype_f'] = rhs_imex_petsc_data  # pass PETSc data type for f
+    description['sweeper_class'] = imex_1st_order  # pass sweeper (see part B)
     description['sweeper_params'] = sweeper_params  # pass sweeper parameters
     description['level_params'] = level_params  # pass level parameters
     description['step_params'] = step_params  # pass step parameters
     description['space_transfer_class'] = mesh_to_mesh_petsc_dmda  # pass spatial transfer class
-    # description['space_transfer_params'] = space_transfer_params  # pass paramters for spatial transfer
+    description['space_transfer_params'] = space_transfer_params  # pass paramters for spatial transfer
 
     # set time parameters
     t0 = 0.0
-    Tend = 1.0
+    Tend = 0.25
 
     # instantiate controller
-    controller = allinclusive_multigrid_MPI(controller_params=controller_params, description=description, comm=time_comm)
-    # controller = allinclusive_multigrid_nonMPI(num_procs=2, controller_params=controller_params, description=description)
+    controller = allinclusive_multigrid_MPI(controller_params=controller_params, description=description,
+                                            comm=time_comm)
 
     # get initial values on finest level
     P = controller.S.levels[0].prob
@@ -117,8 +105,6 @@ def main():
     uex = P.u_exact(Tend)
     err = abs(uex - uend)
 
-    print(err)
-
     # filter statistics by type (number of iterations)
     filtered_stats = filter_stats(stats, type='niter')
 
@@ -143,7 +129,12 @@ def main():
 
     timing = sort_stats(filter_stats(stats, type='timing_run'), sortby='time')
 
-    print(timing)
+    print('Time to solution: %6.4f sec.' % timing[0][1])
+    print('Error vs. PDE solution: %6.4e' % err)
+    print()
+
+    assert err < 2E-04, 'ERROR: did not match error tolerance, got %s' % err
+    assert np.mean(niters) <= 12, 'ERROR: number of iterations is too high, got %s' % np.mean(niters)
 
 
 if __name__ == "__main__":