Skip to content

Commit d1a0d69

Browse files
pancettapancetta
committed
started conv test project
1 parent 5f230ea commit d1a0d69

File tree

6 files changed

+461
-7
lines changed

6 files changed

+461
-7
lines changed

‎projects/AsympConv/PFASST_conv_tests.py‎

Lines changed: 286 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,286 @@
1+
import numpy as np
2+
import pickle
3+
import os
4+
import matplotlib
5+
# matplotlib.use('Agg')
6+
import matplotlib.pyplot as plt
7+
8+
from pySDC.implementations.problem_classes.HeatEquation_1D_FD import heat1d
9+
from pySDC.implementations.problem_classes.AdvectionEquation_1D_FD import advection1d
10+
from pySDC.implementations.datatype_classes.mesh import mesh
11+
from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
12+
from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
13+
from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
14+
from pySDC.implementations.controller_classes.allinclusive_multigrid_nonMPI import allinclusive_multigrid_nonMPI
15+
16+
from pySDC.helpers.stats_helper import filter_stats, sort_stats
17+
18+
19+
def main():
20+
run_diffusion()
21+
run_advection()
22+
23+
24+
def run_diffusion():
25+
"""
26+
A simple test program to test PFASST convergence for the periodic heat equation with random initial data
27+
"""
28+
29+
# initialize level parameters
30+
level_params = dict()
31+
level_params['restol'] = 1E-08
32+
level_params['dt'] = 0.25
33+
level_params['nsweeps'] = [3, 1]
34+
35+
# initialize sweeper parameters
36+
sweeper_params = dict()
37+
sweeper_params['collocation_class'] = CollGaussRadau_Right
38+
sweeper_params['num_nodes'] = [3]
39+
sweeper_params['QI'] = ['LU']
40+
sweeper_params['spread'] = False
41+
42+
# initialize problem parameters
43+
problem_params = dict()
44+
problem_params['freq'] = -1 # frequency for the test value
45+
problem_params['nvars'] = [127, 63] # number of degrees of freedom for each level
46+
47+
# initialize step parameters
48+
step_params = dict()
49+
step_params['maxiter'] = 50
50+
51+
# initialize space transfer parameters
52+
space_transfer_params = dict()
53+
space_transfer_params['rorder'] = 2
54+
space_transfer_params['iorder'] = 2
55+
space_transfer_params['periodic'] = False
56+
57+
# initialize controller parameters
58+
controller_params = dict()
59+
controller_params['logger_level'] = 30
60+
controller_params['predict'] = False
61+
62+
# fill description dictionary for easy step instantiation
63+
description = dict()
64+
description['problem_class'] = heat1d # pass problem class
65+
description['dtype_u'] = mesh # pass data type for u
66+
description['dtype_f'] = mesh # pass data type for f
67+
description['sweeper_class'] = generic_implicit # pass sweeper (see part B)
68+
description['sweeper_params'] = sweeper_params # pass sweeper parameters
69+
description['level_params'] = level_params # pass level parameters
70+
description['step_params'] = step_params # pass step parameters
71+
description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
72+
description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer
73+
74+
# set time parameters
75+
t0 = 0.0
76+
Tend = 4 * level_params['dt']
77+
78+
# set up number of parallel time-steps to run PFASST with
79+
num_proc = 4
80+
81+
results = dict()
82+
83+
for i in range(-3, 12):
84+
ratio = level_params['dt'] / (1.0 / (problem_params['nvars'][0] + 1)) ** 2
85+
86+
problem_params['nu'] = 10.0 ** i / ratio # diffusion coefficient
87+
description['problem_params'] = problem_params # pass problem parameters
88+
89+
out = 'Working on c = %6.4e' % problem_params['nu']
90+
print(out)
91+
cfl = ratio * problem_params['nu']
92+
out = ' CFL number: %4.2e' % cfl
93+
print(out)
94+
95+
# instantiate controller
96+
controller = allinclusive_multigrid_nonMPI(num_procs=num_proc, controller_params=controller_params,
97+
description=description)
98+
99+
# get initial values on finest level
100+
P = controller.MS[0].levels[0].prob
101+
uinit = P.u_exact(t0)
102+
103+
# call main function to get things done...
104+
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
105+
106+
# filter statistics by type (number of iterations)
107+
filtered_stats = filter_stats(stats, type='niter')
108+
109+
# convert filtered statistics to list of iterations count, sorted by process
110+
iter_counts = sort_stats(filtered_stats, sortby='time')
111+
112+
niters = np.array([item[1] for item in iter_counts])
113+
114+
out = ' Mean number of iterations: %4.2f' % np.mean(niters)
115+
print(out)
116+
117+
results[cfl] = np.mean(niters)
118+
119+
file = open('data/results_conv_diffusion.pkl', 'wb')
120+
pickle.dump(results, file)
121+
file.close()
122+
123+
assert os.path.isfile('data/results_conv_diffusion.pkl'), 'ERROR: pickle did not create file'
124+
125+
126+
def run_advection():
127+
"""
128+
A simple test program to test PFASST convergence for the periodic heat equation with random initial data
129+
"""
130+
131+
# initialize level parameters
132+
level_params = dict()
133+
level_params['restol'] = 1E-08
134+
level_params['dt'] = 0.25
135+
level_params['nsweeps'] = [3, 1]
136+
137+
# initialize sweeper parameters
138+
sweeper_params = dict()
139+
sweeper_params['collocation_class'] = CollGaussRadau_Right
140+
sweeper_params['num_nodes'] = [3]
141+
sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
142+
sweeper_params['spread'] = False
143+
144+
# initialize problem parameters
145+
problem_params = dict()
146+
problem_params['freq'] = 64 # frequency for the test value
147+
problem_params['nvars'] = [128, 64] # number of degrees of freedom for each level
148+
problem_params['order'] = 2
149+
problem_params['type'] = 'center'
150+
151+
# initialize step parameters
152+
step_params = dict()
153+
step_params['maxiter'] = 50
154+
155+
# initialize space transfer parameters
156+
space_transfer_params = dict()
157+
space_transfer_params['rorder'] = 2
158+
space_transfer_params['iorder'] = 2
159+
space_transfer_params['periodic'] = True
160+
161+
# initialize controller parameters
162+
controller_params = dict()
163+
controller_params['logger_level'] = 30
164+
controller_params['predict'] = False
165+
166+
# fill description dictionary for easy step instantiation
167+
description = dict()
168+
description['problem_class'] = advection1d # pass problem class
169+
description['dtype_u'] = mesh # pass data type for u
170+
description['dtype_f'] = mesh # pass data type for f
171+
description['sweeper_class'] = generic_implicit # pass sweeper (see part B)
172+
description['sweeper_params'] = sweeper_params # pass sweeper parameters
173+
description['level_params'] = level_params # pass level parameters
174+
description['step_params'] = step_params # pass step parameters
175+
description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
176+
description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer
177+
178+
# set time parameters
179+
t0 = 0.0
180+
Tend = 4 * level_params['dt']
181+
182+
# set up number of parallel time-steps to run PFASST with
183+
num_proc = 4
184+
185+
results = dict()
186+
187+
for i in range(-3, 12):
188+
ratio = level_params['dt'] / (1.0 / (problem_params['nvars'][0] + 1))
189+
190+
problem_params['c'] = 10.0 ** i / ratio # diffusion coefficient
191+
description['problem_params'] = problem_params # pass problem parameters
192+
193+
out = 'Working on nu = %6.4e' % problem_params['c']
194+
print(out)
195+
cfl = ratio * problem_params['c']
196+
out = ' CFL number: %4.2e' % cfl
197+
print(out)
198+
199+
# instantiate controller
200+
controller = allinclusive_multigrid_nonMPI(num_procs=num_proc, controller_params=controller_params,
201+
description=description)
202+
203+
# get initial values on finest level
204+
P = controller.MS[0].levels[0].prob
205+
uinit = P.u_exact(t0)
206+
207+
# call main function to get things done...
208+
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
209+
210+
# filter statistics by type (number of iterations)
211+
filtered_stats = filter_stats(stats, type='niter')
212+
213+
# convert filtered statistics to list of iterations count, sorted by process
214+
iter_counts = sort_stats(filtered_stats, sortby='time')
215+
216+
niters = np.array([item[1] for item in iter_counts])
217+
218+
out = ' Mean number of iterations: %4.2f' % np.mean(niters)
219+
print(out)
220+
221+
results[cfl] = np.mean(niters)
222+
223+
file = open('data/results_conv_advection.pkl', 'wb')
224+
pickle.dump(results, file)
225+
file.close()
226+
227+
assert os.path.isfile('data/results_conv_advection.pkl'), 'ERROR: pickle did not create file'
228+
229+
230+
def plot_results():
231+
232+
file = open('data/results_conv_diffusion.pkl', 'rb')
233+
results_diff = pickle.load(file)
234+
file.close()
235+
236+
file = open('data/results_conv_advection.pkl', 'rb')
237+
results_adv = pickle.load(file)
238+
file.close()
239+
240+
xvalues_diff = sorted(list(results_diff.keys()))
241+
niter_diff = []
242+
for x in xvalues_diff:
243+
niter_diff.append(results_diff[x])
244+
245+
xvalues_adv = sorted(list(results_adv.keys()))
246+
niter_adv = []
247+
for x in xvalues_adv:
248+
niter_adv.append(results_adv[x])
249+
250+
# set up plotting parameters
251+
params = {'legend.fontsize': 20,
252+
'figure.figsize': (12, 8),
253+
'axes.labelsize': 20,
254+
'axes.titlesize': 20,
255+
'xtick.labelsize': 16,
256+
'ytick.labelsize': 16,
257+
'lines.linewidth': 3
258+
}
259+
plt.rcParams.update(params)
260+
261+
# set up figure
262+
plt.figure()
263+
plt.xlabel(r'$\mu$')
264+
plt.ylabel('#iterations')
265+
plt.xlim(min(xvalues_diff + xvalues_adv) / 10, max(xvalues_diff + xvalues_adv) * 10)
266+
plt.ylim(min(niter_diff + niter_adv) - 1, max(niter_diff + niter_adv) + 1)
267+
plt.grid()
268+
269+
# plot
270+
plt.semilogx(xvalues_diff, niter_diff, 'r-', marker='s', markersize=10, label='diffusion')
271+
plt.semilogx(xvalues_adv, niter_adv, 'b-', marker='o', markersize=10, label='advection')
272+
# plt.plot()
273+
274+
plt.legend(loc=1, ncol=1, numpoints=1)
275+
276+
# plt.show()
277+
# save plot, beautify
278+
fname = 'data/conv_test_niter.png'
279+
plt.savefig(fname, rasterized=True, bbox_inches='tight')
280+
281+
assert os.path.isfile(fname), 'ERROR: plotting did not create file'
282+
283+
284+
if __name__ == "__main__":
285+
main()
286+
plot_results()

‎projects/AsympConv/__init__.py‎

Whitespace-only changes.

‎pySDC/core/Sweeper.py‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@ class __Pars(FrozenClass):
3434
def __init__(self, pars):
3535

3636
self.do_coll_update = False
37+
self.spread = True
3738

3839
for k, v in pars.items():
3940
if k is not 'collocation_class':
@@ -138,8 +139,12 @@ def predict(self):
138139

139140
# copy u[0] to all collocation nodes, evaluate RHS
140141
for m in range(1, self.coll.num_nodes + 1):
141-
L.u[m] = P.dtype_u(L.u[0])
142-
L.f[m] = P.eval_f(L.u[m], L.time + L.dt * self.coll.nodes[m - 1])
142+
if self.params.spread:
143+
L.u[m] = P.dtype_u(L.u[0])
144+
L.f[m] = P.eval_f(L.u[m], L.time + L.dt * self.coll.nodes[m - 1])
145+
else:
146+
L.u[m] = P.dtype_u(init=P.init, val=0)
147+
L.f[m] = P.dtype_f(init=P.init, val=0)
143148

144149
# indicate that this level is now ready for sweeps
145150
L.status.unlocked = True

‎pySDC/implementations/problem_classes/AdvectionEquation_1D_FD.py‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,8 @@ def __init__(self, problem_params, dtype_u, dtype_f):
3838
# we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on
3939
if (problem_params['nvars']) % 2 != 0:
4040
raise ProblemError('setup requires nvars = 2^p')
41+
if problem_params['freq'] >= 0 and problem_params['freq'] % 2 != 0:
42+
raise ProblemError('need even number of frequencies due to periodic BCs')
4143

4244
if 'order' not in problem_params:
4345
problem_params['order'] = 1
@@ -175,6 +177,10 @@ def u_exact(self, t):
175177
"""
176178

177179
me = self.dtype_u(self.init)
178-
xvalues = np.array([i * self.dx for i in range(self.params.nvars)])
179-
me.values = np.sin(np.pi * self.params.freq * (xvalues - self.params.c * t))
180+
if self.params.freq >= 0:
181+
xvalues = np.array([i * self.dx for i in range(self.params.nvars)])
182+
me.values = np.sin(np.pi * self.params.freq * (xvalues - self.params.c * t))
183+
else:
184+
np.random.seed(1)
185+
me.values = np.random.rand(self.params.nvars)
180186
return me

0 commit comments

Comments
 (0)