Coverage for pySDC/tutorial/step_7/C_pySDC_with_PETSc.py: 97%
95 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
1import sys
2from pathlib import Path
4import numpy as np
5from mpi4py import MPI
7from pySDC.helpers.stats_helper import get_sorted
9from pySDC.implementations.controller_classes.controller_MPI import controller_MPI
10from pySDC.implementations.problem_classes.HeatEquation_2D_PETSc_forced import heat2d_petsc_forced
11from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
12from pySDC.implementations.transfer_classes.TransferPETScDMDA import mesh_to_mesh_petsc_dmda
15def main():
16 """
17 Program to demonstrate usage of PETSc data structures and spatial parallelization,
18 combined with parallelization in time.
19 """
20 # set MPI communicator
21 comm = MPI.COMM_WORLD
23 world_rank = comm.Get_rank()
24 world_size = comm.Get_size()
26 # split world communicator to create space-communicators
27 if len(sys.argv) >= 2:
28 color = int(world_rank / int(sys.argv[1]))
29 else:
30 color = int(world_rank / 1)
31 space_comm = comm.Split(color=color)
32 space_rank = space_comm.Get_rank()
34 # split world communicator to create time-communicators
35 if len(sys.argv) >= 2:
36 color = int(world_rank % int(sys.argv[1]))
37 else:
38 color = int(world_rank / world_size)
39 time_comm = comm.Split(color=color)
40 time_rank = time_comm.Get_rank()
42 # initialize level parameters
43 level_params = dict()
44 level_params['restol'] = 1e-08
45 level_params['dt'] = 0.125
46 level_params['nsweeps'] = [1]
48 # initialize sweeper parameters
49 sweeper_params = dict()
50 sweeper_params['quad_type'] = 'RADAU-RIGHT'
51 sweeper_params['num_nodes'] = [3]
52 sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
53 sweeper_params['initial_guess'] = 'zero'
55 # initialize problem parameters
56 problem_params = dict()
57 problem_params['nu'] = 1.0 # diffusion coefficient
58 problem_params['freq'] = 2 # frequency for the test value
59 problem_params['cnvars'] = [(65, 65)] # number of degrees of freedom for the coarsest level
60 problem_params['refine'] = [1, 0] # number of refinements
61 problem_params['comm'] = space_comm # pass space-communicator to problem class
62 problem_params['sol_tol'] = 1e-12 # set tolerance to PETSc' linear solver
64 # initialize step parameters
65 step_params = dict()
66 step_params['maxiter'] = 50
68 # initialize space transfer parameters
69 space_transfer_params = dict()
70 space_transfer_params['rorder'] = 2
71 space_transfer_params['iorder'] = 2
72 space_transfer_params['periodic'] = False
74 # initialize controller parameters
75 controller_params = dict()
76 controller_params['logger_level'] = 20 if space_rank == 0 else 99 # set level depending on rank
77 controller_params['dump_setup'] = False
79 # fill description dictionary for easy step instantiation
80 description = dict()
81 description['problem_class'] = heat2d_petsc_forced # pass problem class
82 description['problem_params'] = problem_params # pass problem parameters
83 description['sweeper_class'] = imex_1st_order # pass sweeper (see part B)
84 description['sweeper_params'] = sweeper_params # pass sweeper parameters
85 description['level_params'] = level_params # pass level parameters
86 description['step_params'] = step_params # pass step parameters
87 description['space_transfer_class'] = mesh_to_mesh_petsc_dmda # pass spatial transfer class
88 description['space_transfer_params'] = space_transfer_params # pass parameters for spatial transfer
90 # set time parameters
91 t0 = 0.0
92 Tend = 0.25
94 # instantiate controller
95 controller = controller_MPI(controller_params=controller_params, description=description, comm=time_comm)
97 # get initial values on finest level
98 P = controller.S.levels[0].prob
99 uinit = P.u_exact(t0)
101 # call main function to get things done...
102 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
104 # compute exact solution and compare
105 uex = P.u_exact(Tend)
106 err = abs(uex - uend)
108 # filter statistics by type (number of iterations)
109 iter_counts = get_sorted(stats, type='niter', sortby='time')
111 niters = np.array([item[1] for item in iter_counts])
113 # limit output to space-rank 0 (as before when setting the logger level)
114 if space_rank == 0:
115 if len(sys.argv) == 3:
116 fname = str(sys.argv[2])
117 else:
118 fname = 'step_7_C_out.txt'
119 Path("data").mkdir(parents=True, exist_ok=True)
120 f = open('data/' + fname, 'a+')
122 out = 'This is time-rank %i...' % time_rank
123 f.write(out + '\n')
124 print(out)
126 # compute and print statistics
127 for item in iter_counts:
128 out = 'Number of iterations for time %4.2f: %2i' % item
129 f.write(out + '\n')
130 print(out)
132 out = ' Mean number of iterations: %4.2f' % np.mean(niters)
133 f.write(out + '\n')
134 print(out)
135 out = ' Range of values for number of iterations: %2i ' % np.ptp(niters)
136 f.write(out + '\n')
137 print(out)
138 out = ' Position of max/min number of iterations: %2i -- %2i' % (
139 int(np.argmax(niters)),
140 int(np.argmin(niters)),
141 )
142 f.write(out + '\n')
143 print(out)
144 out = ' Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters)))
145 f.write(out + '\n')
146 print(out)
148 timing = get_sorted(stats, type='timing_run', sortby='time')
150 out = 'Time to solution: %6.4f sec.' % timing[0][1]
151 f.write(out + '\n')
152 print(out)
153 out = 'Error vs. PDE solution: %6.4e' % err
154 f.write(out + '\n')
155 print(out)
157 f.close()
159 assert err < 2e-04, 'ERROR: did not match error tolerance, got %s' % err
160 assert np.mean(niters) <= 12, 'ERROR: number of iterations is too high, got %s' % np.mean(niters)
162 space_comm.Free()
163 time_comm.Free()
166if __name__ == "__main__":
167 main()