Coverage for pySDC/projects/matrixPFASST/compare_to_propagator.py: 100%
135 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
1import numpy as np
2from pathlib import Path
4from pySDC.helpers.stats_helper import get_sorted
6from pySDC.implementations.problem_classes.AdvectionEquation_ND_FD import advectionNd
7from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
8from pySDC.implementations.problem_classes.TestEquation_0D import testequation0d
9from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
10from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
11from pySDC.implementations.transfer_classes.TransferMesh_NoCoarse import mesh_to_mesh as mesh_to_mesh_nocoarse
12from pySDC.projects.matrixPFASST.controller_matrix_nonMPI import controller_matrix_nonMPI
15def diffusion_setup(par=0.0):
16 """
17 Setup routine for advection test
19 Args:
20 par (float): parameter for controlling stiffness
21 """
22 # initialize level parameters
23 level_params = dict()
24 level_params['restol'] = 1e-08
25 level_params['dt'] = 0.25
26 level_params['nsweeps'] = [3, 1]
28 # initialize sweeper parameters
29 sweeper_params = dict()
30 sweeper_params['quad_type'] = 'RADAU-RIGHT'
31 sweeper_params['num_nodes'] = 3
32 sweeper_params['QI'] = 'LU'
33 sweeper_params['initial_guess'] = 'spread'
35 # initialize problem parameters
36 problem_params = dict()
37 problem_params['nu'] = par # diffusion coefficient
38 problem_params['freq'] = 4 # frequency for the test value
39 problem_params['nvars'] = [127] # number of degrees of freedom for each level
40 problem_params['bc'] = 'dirichlet-zero' # boundary conditions
42 # initialize step parameters
43 step_params = dict()
44 step_params['maxiter'] = 50
46 # initialize space transfer parameters
47 space_transfer_params = dict()
48 space_transfer_params['rorder'] = 2
49 space_transfer_params['iorder'] = 2
51 # initialize controller parameters
52 controller_params = dict()
53 controller_params['logger_level'] = 30
55 # fill description dictionary for easy step instantiation
56 description = dict()
57 description['problem_class'] = heatNd_unforced # pass problem class
58 description['problem_params'] = problem_params # pass problem parameters
59 description['sweeper_class'] = generic_implicit # pass sweeper
60 description['sweeper_params'] = sweeper_params # pass sweeper parameters
61 description['level_params'] = level_params # pass level parameters
62 description['step_params'] = step_params # pass step parameters
63 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
64 description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer
66 return description, controller_params
69def advection_setup(par=0.0):
70 """
71 Setup routine for advection test
73 Args:
74 par (float): parameter for controlling stiffness
75 """
76 # initialize level parameters
77 level_params = dict()
78 level_params['restol'] = 1e-08
79 level_params['dt'] = 0.25
80 level_params['nsweeps'] = [3, 1]
82 # initialize sweeper parameters
83 sweeper_params = dict()
84 sweeper_params['quad_type'] = 'RADAU-RIGHT'
85 sweeper_params['num_nodes'] = [3]
86 sweeper_params['QI'] = ['LU']
87 sweeper_params['initial_guess'] = 'spread'
89 # initialize problem parameters
90 problem_params = dict()
91 problem_params['c'] = par
92 problem_params['freq'] = 4 # frequency for the test value
93 problem_params['nvars'] = [128, 64] # number of degrees of freedom for each level
94 problem_params['order'] = 2
95 problem_params['stencil_type'] = 'center'
96 problem_params['bc'] = 'periodic' # boundary conditions
98 # initialize step parameters
99 step_params = dict()
100 step_params['maxiter'] = 50
102 # initialize space transfer parameters
103 space_transfer_params = dict()
104 space_transfer_params['rorder'] = 2
105 space_transfer_params['iorder'] = 2
106 space_transfer_params['periodic'] = True
108 # initialize controller parameters
109 controller_params = dict()
110 controller_params['logger_level'] = 30
112 # fill description dictionary for easy step instantiation
113 description = dict()
114 description['problem_class'] = advectionNd # pass problem class
115 description['problem_params'] = problem_params
116 description['sweeper_class'] = generic_implicit # pass sweeper (see part B)
117 description['sweeper_params'] = sweeper_params # pass sweeper parameters
118 description['level_params'] = level_params # pass level parameters
119 description['step_params'] = step_params # pass step parameters
120 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
121 description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer
123 return description, controller_params
126def scalar_equation_setup():
127 """
128 Setup routine for the test equation
130 Args:
131 par (float): parameter for controlling stiffness
132 """
133 # initialize level parameters
134 level_params = dict()
135 level_params['restol'] = 1e-08
136 level_params['dt'] = 0.25
137 level_params['nsweeps'] = [3, 1]
139 # initialize sweeper parameters
140 sweeper_params = dict()
141 sweeper_params['quad_type'] = 'RADAU-RIGHT'
142 sweeper_params['num_nodes'] = [3, 2]
143 sweeper_params['QI'] = 'LU'
144 sweeper_params['initial_guess'] = 'spread'
146 # initialize problem parameters
147 problem_params = dict()
148 problem_params['u0'] = 1.0 # initial value (for all instances)
149 # use single values like this...
150 # problem_params['lambdas'] = [[-1.0]]
151 # .. or a list of values like this ...
152 # problem_params['lambdas'] = [[-1.0, -2.0, 1j, -1j]]
153 # .. or a whole block of values like this
154 ilim_left = -11
155 ilim_right = 0
156 rlim_left = 0
157 rlim_right = 11
158 ilam = 1j * np.logspace(ilim_left, ilim_right, 11)
159 rlam = -1 * np.logspace(rlim_left, rlim_right, 11)
160 lambdas = []
161 for rl in rlam:
162 for il in ilam:
163 lambdas.append(rl + il)
164 problem_params['lambdas'] = [lambdas]
165 # note: PFASST will do all of those at once, but without interaction (realized via diagonal matrix).
166 # The propagation matrix will be diagonal too, corresponding to the respective lambda value.
168 # initialize step parameters
169 step_params = dict()
170 step_params['maxiter'] = 50
172 # initialize controller parameters
173 controller_params = dict()
174 controller_params['logger_level'] = 30
176 # fill description dictionary for easy step instantiation
177 description = dict()
178 description['problem_class'] = testequation0d # pass problem class
179 description['problem_params'] = problem_params # pass problem parameters
180 description['sweeper_class'] = generic_implicit # pass sweeper
181 description['sweeper_params'] = sweeper_params # pass sweeper parameters
182 description['level_params'] = level_params # pass level parameters
183 description['step_params'] = step_params # pass step parameters
184 description['space_transfer_class'] = mesh_to_mesh_nocoarse # pass spatial transfer class
185 description['space_transfer_params'] = dict() # pass paramters for spatial transfer
187 return description, controller_params
190def compare_controllers(type=None, par=0.0, f=None):
191 """
192 A simple test program to compare PFASST runs with matrix-based and matrix-free controllers
194 Args:
195 type (str): setup type
196 par (float): parameter for controlling stiffness
197 f: file handler
198 """
200 # set time parameters
201 t0 = 0.0
202 Tend = 1.0
204 if type == 'diffusion':
205 description, controller_params = diffusion_setup(par)
206 elif type == 'advection':
207 description, controller_params = advection_setup(par)
208 elif type == 'testequation':
209 description, controller_params = scalar_equation_setup()
210 else:
211 raise ValueError('No valis setup type provided, aborting..')
213 out = '\nWorking with %s setup and parameter %3.1e..' % (type, par)
214 f.write(out + '\n')
215 print(out)
217 # instantiate controller
218 controller = controller_matrix_nonMPI(num_procs=4, controller_params=controller_params, description=description)
219 # get initial values on finest level
220 P = controller.MS[0].levels[0].prob
221 uinit = P.u_exact(t0)
222 uex = P.u_exact(Tend)
224 # this is where the iteration is happening
225 uend_mat, stats_mat = controller.run(u0=uinit, t0=t0, Tend=Tend)
227 # filter statistics by type (number of iterations)
228 iter_counts_mat = get_sorted(stats_mat, type='niter', sortby='time')
230 out = ' Iteration counts for matrix-based version: %s' % iter_counts_mat
231 f.write(out + '\n')
232 print(out)
234 # filter only iteration counts and check for equality
235 niters = [item[1] for item in iter_counts_mat]
236 assert niters.count(niters[0]) == len(niters), 'ERROR: not all time-steps have the same number of iterations'
237 niter = niters[0]
239 # build propagation matrix using the prescribed number of iterations (or any other, if needed)
240 prop = controller.build_propagation_matrix(niter=niter)
242 err_prop_ex = np.linalg.norm(prop.dot(uinit) - uex)
243 err_mat_ex = np.linalg.norm(uend_mat - uex)
244 out = ' Error (mat/prop) vs. exact solution: %6.4e -- %6.4e' % (err_mat_ex, err_prop_ex)
245 f.write(out + '\n')
246 print(out)
247 err_mat_prop = np.linalg.norm(prop.dot(uinit) - uend_mat)
248 out = ' Difference between matrix-PFASST and propagator: %6.4e' % err_mat_prop
249 f.write(out + '\n')
250 print(out)
252 assert err_mat_prop < 2.0e-14, (
253 'ERROR: difference between matrix-based and propagator result is too large, got %s' % err_mat_prop
254 )
257def main():
258 par_list = [1e-02, 1.0, 1e02]
260 Path("data").mkdir(parents=True, exist_ok=True)
261 f = open('data/comparison_matrix_vs_propagator_detail.txt', 'w')
262 for par in par_list:
263 compare_controllers(type='diffusion', par=par, f=f)
264 compare_controllers(type='advection', par=par, f=f)
265 compare_controllers(type='testequation', par=0.0, f=f)
266 f.close()
269if __name__ == "__main__":
270 main()