Skip to content

Commit d22fec0

Browse files
pancettapancetta
committed
Merge pull request #84 from danielru/add/stab_vs_k
Add/stab vs k
2 parents 9d2d65a + 3a8fee9 commit d22fec0

File tree

4 files changed

+145
-45
lines changed

4 files changed

+145
-45
lines changed

‎examples/fwsw/plot_stab_vs_k.py‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
from pySDC import Level as lvl
2+
from pySDC import Hooks as hookclass
3+
from pySDC import CollocationClasses as collclass
4+
from pySDC import Step as stepclass
5+
6+
from pySDC.datatype_classes.complex_mesh import mesh, rhs_imex_mesh
7+
from pySDC.sweeper_classes.imex_1st_order import imex_1st_order as imex
8+
from examples.fwsw.ProblemClass import swfw_scalar
9+
import numpy as np
10+
11+
from pylab import rcParams
12+
import matplotlib.pyplot as plt
13+
from matplotlib.patches import Polygon
14+
from subprocess import call
15+
from matplotlib.ticker import ScalarFormatter
16+
17+
if __name__ == "__main__":
18+
mvals = [2, 3, 4]
19+
kvals = np.arange(1,10)
20+
lambda_fast = 10j
21+
22+
### PLOT EITHER FOR lambda_slow = 1 (resolved) OR lambda_slow = 4 (unresolved)
23+
slow_resolved = False
24+
if slow_resolved:
25+
lambda_slow = 1j
26+
else:
27+
lambda_slow = 4j
28+
stabval = np.zeros((np.size(mvals), np.size(kvals)))
29+
30+
for i in range(0,np.size(mvals)):
31+
pparams = {}
32+
# the following are not used in the computation
33+
pparams['lambda_s'] = np.array([0.0])
34+
pparams['lambda_f'] = np.array([0.0])
35+
pparams['u0'] = 1.0
36+
swparams = {}
37+
# swparams['collocation_class'] = collclass.CollGaussLobatto
38+
# swparams['collocation_class'] = collclass.CollGaussLegendre
39+
swparams['collocation_class'] = collclass.CollGaussRadau_Right
40+
swparams['num_nodes'] = mvals[i]
41+
do_coll_update = True
42+
43+
#
44+
# ...this is functionality copied from test_imexsweeper. Ideally, it should be available in one place.
45+
#
46+
step = stepclass.step(params={})
47+
L = lvl.level(problem_class=swfw_scalar, problem_params=pparams, dtype_u=mesh, dtype_f=rhs_imex_mesh, sweeper_class=imex, sweeper_params=swparams, level_params={}, hook_class=hookclass.hooks, id="stability")
48+
step.register_level(L)
49+
step.status.dt = 1.0 # Needs to be 1.0, change dt through lambdas
50+
step.status.time = 0.0
51+
u0 = step.levels[0].prob.u_exact(step.status.time)
52+
step.init_step(u0)
53+
nnodes = step.levels[0].sweep.coll.num_nodes
54+
level = step.levels[0]
55+
problem = level.prob
56+
57+
QE = level.sweep.QE[1:,1:]
58+
QI = level.sweep.QI[1:,1:]
59+
Q = level.sweep.coll.Qmat[1:,1:]
60+
61+
LHS, RHS = level.sweep.get_scalar_problems_sweeper_mats( lambdas = [ lambda_fast, lambda_slow ] )
62+
63+
for k in range(0, np.size(kvals)):
64+
Kmax = kvals[k]
65+
Mat_sweep = level.sweep.get_scalar_problems_manysweep_mat( nsweeps = Kmax, lambdas = [ lambda_fast, lambda_slow ] )
66+
if do_coll_update:
67+
stab_fh = 1.0 + (lambda_fast + lambda_slow)*level.sweep.coll.weights.dot(Mat_sweep.dot(np.ones(nnodes)))
68+
else:
69+
q = np.zeros(nnodes)
70+
q[nnodes-1] = 1.0
71+
stab_fh = q.dot(Mat_sweep.dot(np.ones(nnodes)))
72+
stabval[i,k] = np.absolute(stab_fh)
73+
74+
rcParams['figure.figsize'] = 2.5, 2.5
75+
fig = plt.figure()
76+
fs = 8
77+
plt.plot(kvals, stabval[0,:], 'o-', color='b', label=("M=%2i" % mvals[0]), markersize=fs-2)
78+
plt.plot(kvals, stabval[1,:], 's-', color='r', label=("M=%2i" % mvals[1]), markersize=fs-2)
79+
plt.plot(kvals, stabval[2,:], 'd-', color='g', label=("M=%2i" % mvals[2]), markersize=fs-2)
80+
plt.plot(kvals, 1.0+0.0*kvals, '--', color='k')
81+
plt.xlabel('Number of iterations K', fontsize=fs)
82+
plt.ylabel(r'Modulus of stability function $\left| R \right|$', fontsize=fs)
83+
plt.ylim([0.0, 1.2])
84+
if slow_resolved:
85+
plt.legend(loc='upper right', fontsize=fs, prop={'size':fs})
86+
else:
87+
plt.legend(loc='lower left', fontsize=fs, prop={'size':fs})
88+
89+
plt.gca().get_xaxis().get_major_formatter().labelOnlyBase = False
90+
plt.gca().get_xaxis().set_major_formatter(ScalarFormatter())
91+
#plt.show()
92+
if slow_resolved:
93+
filename = 'stab_vs_k_resolved.pdf'
94+
else:
95+
filename = 'stab_vs_k_unresolved.pdf'
96+
97+
fig.savefig(filename, bbox_inches='tight')
98+
call(["pdfcrop", filename, filename])

‎examples/fwsw/plot_stablimits.py‎ renamed to ‎examples/fwsw/plot_stab_vs_m.py‎

Lines changed: 34 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -12,33 +12,25 @@
1212
import matplotlib.pyplot as plt
1313
from matplotlib.patches import Polygon
1414
from subprocess import call
15-
16-
def get_stab_function(LHS, RHS, Kmax, lambd, do_coll_update):
17-
Pinv = np.linalg.inv(LHS)
18-
Mat_sweep = np.linalg.matrix_power(Pinv.dot(RHS), Kmax)
19-
for l in range(0,Kmax):
20-
Mat_sweep = Mat_sweep + np.linalg.matrix_power(Pinv.dot(RHS),k).dot(Pinv)
21-
if do_coll_update:
22-
stab_fh = 1.0 + lambd*level.sweep.coll.weights.dot(Mat_sweep.dot(np.ones(nnodes)))
23-
else:
24-
q = np.zeros(nnodes)
25-
q[nnodes-1] = 1.0
26-
stab_fh = q.dot(Mat_sweep.dot(np.ones(nnodes)))
27-
return stab_fh
15+
from matplotlib.ticker import ScalarFormatter
2816

2917
if __name__ == "__main__":
30-
mvals = [4]
31-
kvals = np.arange(2,9)
32-
lambdaratio = [1, 5, 50]
33-
stabval = np.zeros((np.size(mvals), np.size(lambdaratio), np.size(kvals)))
18+
mvals = np.arange(2,10)
19+
kvals = [3, 5, 7]
20+
lambda_fast = 10j
21+
lambda_slow = 3j
22+
stabval = np.zeros((np.size(kvals), np.size(mvals)))
3423

3524
for i in range(0,np.size(mvals)):
3625
pparams = {}
3726
# the following are not used in the computation
3827
pparams['lambda_s'] = np.array([0.0])
3928
pparams['lambda_f'] = np.array([0.0])
4029
pparams['u0'] = 1.0
30+
4131
swparams = {}
32+
# swparams['collocation_class'] = collclass.CollGaussLobatto
33+
# swparams['collocation_class'] = collclass.CollGaussLegendre
4234
swparams['collocation_class'] = collclass.CollGaussRadau_Right
4335
swparams['num_nodes'] = mvals[i]
4436
do_coll_update = True
@@ -60,28 +52,34 @@ def get_stab_function(LHS, RHS, Kmax, lambd, do_coll_update):
6052
QE = level.sweep.QE[1:,1:]
6153
QI = level.sweep.QI[1:,1:]
6254
Q = level.sweep.coll.Qmat[1:,1:]
55+
LHS, RHS = level.sweep.get_scalar_problems_sweeper_mats( lambdas = [ lambda_fast, lambda_slow ] )
6356

64-
for j in range(0,np.size(lambdaratio)):
65-
lambda_slow = 0.1*1j
66-
lambda_fast = lambdaratio[j]*lambda_slow
67-
68-
LHS = np.eye(nnodes) - step.status.dt*( lambda_fast*QI + lambda_slow*QE )
69-
RHS = step.status.dt*( (lambda_fast+lambda_slow)*Q - (lambda_fast*QI + lambda_slow*QE) )
70-
71-
for k in range(0, np.size(kvals)):
72-
Kmax = kvals[k]
73-
stab_fh = get_stab_function(LHS, RHS, kvals[k], lambda_fast+lambda_slow, do_coll_update)
74-
stabval[i,j,k] = np.absolute(stab_fh)
57+
for k in range(0,np.size(kvals)):
58+
Kmax = kvals[k]
59+
Mat_sweep = level.sweep.get_scalar_problems_manysweep_mat( nsweeps = Kmax, lambdas = [ lambda_fast, lambda_slow ] )
60+
if do_coll_update:
61+
stab_fh = 1.0 + (lambda_fast + lambda_slow)*level.sweep.coll.weights.dot(Mat_sweep.dot(np.ones(nnodes)))
62+
else:
63+
q = np.zeros(nnodes)
64+
q[nnodes-1] = 1.0
65+
stab_fh = q.dot(Mat_sweep.dot(np.ones(nnodes)))
66+
stabval[k,i] = np.absolute(stab_fh)
7567

68+
rcParams['figure.figsize'] = 2.5, 2.5
7669
fig = plt.figure()
7770
fs = 8
78-
plt.plot(kvals, stabval[0,0,:], 'o-', color='b', label=("Ratio: %3.0f" % lambdaratio[0]))
79-
plt.plot(kvals, stabval[0,1,:], 's-', color='r', label=("Ratio: %3.0f" % lambdaratio[1]))
80-
plt.plot(kvals, stabval[0,2,:], 'd-', color='g', label=("Ratio: %3.0f" % lambdaratio[2]))
81-
plt.plot(kvals, 1.0+0.0*kvals, '--', color='k')
71+
plt.plot(mvals, stabval[0,:], 'o-', color='b', label=(r"K=%1i" % kvals[0]))
72+
plt.plot(mvals, stabval[1,:], 's-', color='r', label=(r"K=%1i" % kvals[1]))
73+
plt.plot(mvals, stabval[2,:], 'd-', color='g', label=(r"K=%1i" % kvals[2]))
74+
plt.plot(mvals, 1.0+0.0*mvals, '--', color='k')
75+
plt.xlabel('Number of nodes M', fontsize=fs)
76+
plt.ylabel(r'Modulus of stability function $\left| R \right|$', fontsize=fs)
77+
plt.ylim([0.0, 1.8])
8278
plt.legend(loc='lower right', fontsize=fs, prop={'size':fs})
83-
84-
#plt.plot(kvals, stabval[1,0,:], '-', color='r')
85-
#plt.plot(kvals, stabval[1,1,:], '--', color='r')
86-
#plt.plot(kvals, stabval[1,2,:], '-.', color='r')
79+
plt.gca().get_xaxis().get_major_formatter().labelOnlyBase = False
80+
plt.gca().get_xaxis().set_major_formatter(ScalarFormatter())
8781
plt.show()
82+
83+
# filename = 'stablimit-M'+str(mvals[0])+'.pdf'
84+
# fig.savefig(filename, bbox_inches='tight')
85+
# call(["pdfcrop", filename, filename])

‎examples/fwsw/plot_stability.py‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -31,14 +31,14 @@
3131
swparams = {}
3232
### SET TYPE OF QUADRATURE NODES ###
3333
#swparams['collocation_class'] = collclass.CollGaussLobatto
34-
swparams['collocation_class'] = collclass.CollGaussLegendre
35-
#swparams['collocation_class'] = collclass.CollGaussRadau_Right
34+
#swparams['collocation_class'] = collclass.CollGaussLegendre
35+
swparams['collocation_class'] = collclass.CollGaussRadau_Right
3636

3737
### SET NUMBER OF QUADRATURE NODES ###
3838
swparams['num_nodes'] = 3
3939

4040
### SET NUMBER OF ITERATIONS - SET K=0 FOR COLLOCATION SOLUTION ###
41-
K = 4
41+
K = 5
4242

4343
do_coll_update = True
4444

@@ -89,6 +89,9 @@
8989
# levels = np.array([1.0])
9090
CS1 = plt.contour(lambda_s.imag, lambda_f.imag, np.absolute(stab), levels, colors='k', linestyles='dashed')
9191
CS2 = plt.contour(lambda_s.imag, lambda_f.imag, np.absolute(stab), [1.0], colors='k')
92+
# Set markers at points used in plot_stab_vs_k
93+
plt.plot(4, 10, 'x', color='k', markersize=fs-4)
94+
plt.plot(1, 10, 'x', color='k', markersize=fs-4)
9295
plt.clabel(CS1, inline=True, fmt='%3.2f', fontsize=fs-2)
9396
manual_locations = [(1.5, 2.5)]
9497
if K>0: # for K=0 and no 1.0 isoline, this crashes Matplotlib for somer reason

‎examples/ruprecht_speck_2016.README‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,10 @@
11
Figure 1: fwsw/plot_stifflimit_specrad.py
22
Figure 2: fwsw/plot_stability.py
3-
Figure 3: acoustic_1d_imex/plot_dispersion.py
4-
Figure 4 (left): acoustic_1d_imex/runconvergence.py and then plotconvergence.py
5-
Figure 4 (right): acoustic_1d_imex/runitererror.py
6-
Figure 5: acoustic_1d_imex/runmultiscale.py
7-
Figure 6: boussinesq_2d_imex/rungmrescounter.py and then plotgmrescounter.py for sparams['maxiter']=4 and dirk_order=4
3+
Figure 3: fwsw/plot_stab_vs_k.py
4+
Figure 4: acoustic_1d_imex/plot_dispersion.py
5+
Figure 5 (left): acoustic_1d_imex/runconvergence.py and then plotconvergence.py
6+
Figure 5 (right): acoustic_1d_imex/runitererror.py
7+
Figure 6: acoustic_1d_imex/runmultiscale.py
8+
Figure 7: boussinesq_2d_imex/rungmrescounter.py and then plotgmrescounter.py for sparams['maxiter']=4 and dirk_order=4
89

9-
Table 1: boussinesq_2d_imex/rungmrescounter.py and then plotgmrescounter.py
10+
Table 1: boussinesq_2d_imex/rungmrescounter.py and then plotgmrescounter.py

0 commit comments

Comments
 (0)