Coverage for pySDC/projects/matrixPFASST/compare_to_matrixbased.py: 100%
142 statements
« prev ^ index » next coverage.py v7.6.7, created at 2024-11-16 14:51 +0000
1import numpy as np
2from pathlib import Path
4from pySDC.helpers.stats_helper import get_sorted
6from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
7from pySDC.implementations.problem_classes.AdvectionEquation_ND_FD import advectionNd
8from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
9from pySDC.implementations.problem_classes.TestEquation_0D import testequation0d
10from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
11from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
12from pySDC.implementations.transfer_classes.TransferMesh_NoCoarse import mesh_to_mesh as mesh_to_mesh_nocoarse
13from pySDC.projects.matrixPFASST.controller_matrix_nonMPI import controller_matrix_nonMPI
16def diffusion_setup(par=0.0):
17 """
18 Setup routine for advection test
20 Args:
21 par (float): parameter for controlling stiffness
22 """
23 # initialize level parameters
24 level_params = dict()
25 level_params['restol'] = 1e-08
26 level_params['dt'] = 0.25
27 level_params['nsweeps'] = [3, 1]
29 # initialize sweeper parameters
30 sweeper_params = dict()
31 sweeper_params['quad_type'] = 'RADAU-RIGHT'
32 sweeper_params['num_nodes'] = 3
33 sweeper_params['QI'] = 'LU' # For the IMEX sweeper, the LU-trick can be activated for the implicit part
34 sweeper_params['initial_guess'] = 'spread'
36 # initialize problem parameters
37 problem_params = dict()
38 problem_params['nu'] = par # diffusion coefficient
39 problem_params['freq'] = 4 # frequency for the test value
40 problem_params['nvars'] = [127, 63] # number of degrees of freedom for each level
41 problem_params['bc'] = 'dirichlet-zero' # boundary conditions
43 # initialize step parameters
44 step_params = dict()
45 step_params['maxiter'] = 50
47 # initialize space transfer parameters
48 space_transfer_params = dict()
49 space_transfer_params['rorder'] = 2
50 space_transfer_params['iorder'] = 2
52 # initialize controller parameters
53 controller_params = dict()
54 controller_params['logger_level'] = 30
55 controller_params['all_to_done'] = True
57 # fill description dictionary for easy step instantiation
58 description = dict()
59 description['problem_class'] = heatNd_unforced # pass problem class
60 description['problem_params'] = problem_params # pass problem parameters
61 description['sweeper_class'] = generic_implicit # pass sweeper
62 description['sweeper_params'] = sweeper_params # pass sweeper parameters
63 description['level_params'] = level_params # pass level parameters
64 description['step_params'] = step_params # pass step parameters
65 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
66 description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer
68 return description, controller_params
71def advection_setup(par=0.0):
72 """
73 Setup routine for advection test
75 Args:
76 par (float): parameter for controlling stiffness
77 """
78 # initialize level parameters
79 level_params = dict()
80 level_params['restol'] = 1e-08
81 level_params['dt'] = 0.25
82 level_params['nsweeps'] = [3, 1]
84 # initialize sweeper parameters
85 sweeper_params = dict()
86 sweeper_params['quad_type'] = 'RADAU-RIGHT'
87 sweeper_params['num_nodes'] = [3]
88 sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
89 sweeper_params['initial_guess'] = 'spread'
91 # initialize problem parameters
92 problem_params = dict()
93 problem_params['c'] = par
94 problem_params['freq'] = 4 # frequency for the test value
95 problem_params['nvars'] = [128, 64] # number of degrees of freedom for each level
96 problem_params['order'] = 2
97 problem_params['stencil_type'] = 'center'
98 problem_params['bc'] = 'periodic' # boundary conditions
100 # initialize step parameters
101 step_params = dict()
102 step_params['maxiter'] = 50
104 # initialize space transfer parameters
105 space_transfer_params = dict()
106 space_transfer_params['rorder'] = 2
107 space_transfer_params['iorder'] = 2
108 space_transfer_params['periodic'] = True
110 # initialize controller parameters
111 controller_params = dict()
112 controller_params['logger_level'] = 30
113 controller_params['all_to_done'] = True
115 # fill description dictionary for easy step instantiation
116 description = dict()
117 description['problem_class'] = advectionNd # pass problem class
118 description['problem_params'] = problem_params
119 description['sweeper_class'] = generic_implicit # pass sweeper (see part B)
120 description['sweeper_params'] = sweeper_params # pass sweeper parameters
121 description['level_params'] = level_params # pass level parameters
122 description['step_params'] = step_params # pass step parameters
123 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
124 description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer
126 return description, controller_params
129def testequation_setup():
130 """
131 Setup routine for the test equation
133 Args:
134 par (float): parameter for controlling stiffness
135 """
136 # initialize level parameters
137 level_params = dict()
138 level_params['restol'] = 1e-08
139 level_params['dt'] = 0.25
140 level_params['nsweeps'] = [3, 1]
142 # initialize sweeper parameters
143 sweeper_params = dict()
144 sweeper_params['quad_type'] = 'RADAU-RIGHT'
145 sweeper_params['num_nodes'] = [3, 2]
146 sweeper_params['QI'] = 'LU'
147 sweeper_params['initial_guess'] = 'spread'
149 # initialize problem parameters
150 problem_params = dict()
151 problem_params['u0'] = 1.0 # initial value (for all instances)
152 # use single values like this...
153 # problem_params['lambdas'] = [[-1.0]]
154 # .. or a list of values like this ...
155 # problem_params['lambdas'] = [[-1.0, -2.0, 1j, -1j]]
156 # .. or a whole block of values like this
157 ilim_left = -11
158 ilim_right = 0
159 rlim_left = 0
160 rlim_right = 11
161 ilam = 1j * np.logspace(ilim_left, ilim_right, 11)
162 rlam = -1 * np.logspace(rlim_left, rlim_right, 11)
163 lambdas = []
164 for rl in rlam:
165 for il in ilam:
166 lambdas.append(rl + il)
167 problem_params['lambdas'] = [lambdas]
168 # note: PFASST will do all of those at once, but without interaction (realized via diagonal matrix).
169 # The propagation matrix will be diagonal too, corresponding to the respective lambda value.
171 # initialize step parameters
172 step_params = dict()
173 step_params['maxiter'] = 50
175 # initialize controller parameters
176 controller_params = dict()
177 controller_params['logger_level'] = 30
178 controller_params['all_to_done'] = True
180 # fill description dictionary for easy step instantiation
181 description = dict()
182 description['problem_class'] = testequation0d # pass problem class
183 description['problem_params'] = problem_params # pass problem parameters
184 description['sweeper_class'] = generic_implicit # pass sweeper
185 description['sweeper_params'] = sweeper_params # pass sweeper parameters
186 description['level_params'] = level_params # pass level parameters
187 description['step_params'] = step_params # pass step parameters
188 description['space_transfer_class'] = mesh_to_mesh_nocoarse # pass spatial transfer class
189 description['space_transfer_params'] = dict() # pass paramters for spatial transfer
191 return description, controller_params
194def compare_controllers(type=None, par=0.0, f=None):
195 """
196 A simple test program to compare PFASST runs with matrix-based and matrix-free controllers
198 Args:
199 type (str): setup type
200 par (float) parameter for controlling stiffness
201 f: file handler
202 """
204 # set time parameters
205 t0 = 0.0
206 Tend = 1.0
208 if type == 'diffusion':
209 description, controller_params = diffusion_setup(par)
210 elif type == 'advection':
211 description, controller_params = advection_setup(par)
212 elif type == 'testequation':
213 description, controller_params = testequation_setup()
214 else:
215 raise ValueError('No valis setup type provided, aborting..')
217 out = '\nWorking with %s setup and parameter %3.1e..' % (type, par)
218 f.write(out + '\n')
219 print(out)
221 # instantiate controller
222 controller_mat = controller_matrix_nonMPI(num_procs=4, controller_params=controller_params, description=description)
224 controller_nomat = controller_nonMPI(num_procs=4, controller_params=controller_params, description=description)
226 # get initial values on finest level
227 P = controller_nomat.MS[0].levels[0].prob
228 uinit = P.u_exact(t0)
229 uex = P.u_exact(Tend)
231 # this is where the iteration is happening
232 uend_mat, stats_mat = controller_mat.run(u0=uinit, t0=t0, Tend=Tend)
233 uend_nomat, stats_nomat = controller_nomat.run(u0=uinit, t0=t0, Tend=Tend)
235 diff = abs(uend_mat - uend_nomat)
237 err_mat = abs(uend_mat - uex)
238 err_nomat = abs(uend_nomat - uex)
240 out = ' Error (mat/nomat) vs. exact solution: %6.4e -- %6.4e' % (err_mat, err_nomat)
241 f.write(out + '\n')
242 print(out)
243 out = ' Difference between both results: %6.4e' % diff
244 f.write(out + '\n')
245 print(out)
247 assert diff < 2.3e-15, 'ERROR: difference between matrix-based and matrix-free result is too large, got %s' % diff
249 # get and convert statistics to list of iterations count, sorted by process
250 iter_counts_mat = get_sorted(stats_mat, type='niter', sortby='time')
251 iter_counts_nomat = get_sorted(stats_nomat, type='niter', sortby='time')
253 out = ' Iteration counts for matrix-based version: %s' % iter_counts_mat
254 f.write(out + '\n')
255 print(out)
256 out = ' Iteration counts for matrix-free version: %s' % iter_counts_nomat
257 f.write(out + '\n')
258 print(out)
260 assert (
261 iter_counts_nomat == iter_counts_mat
262 ), 'ERROR: number of iterations differ between matrix-based and matrix-free controller'
265def main():
266 par_list = [1e-02, 1.0, 1e02]
268 Path("data").mkdir(parents=True, exist_ok=True)
269 f = open('data/comparison_matrix_vs_nomat_detail.txt', 'w')
270 for par in par_list:
271 compare_controllers(type='diffusion', par=par, f=f)
272 compare_controllers(type='advection', par=par, f=f)
273 compare_controllers(type='testequation', par=0.0, f=f)
274 f.close()
277if __name__ == "__main__":
278 main()