Skip to content

Commit bc1cb58

Browse files
pancettapancetta
committed
Merge pull request #69 from danielru/add/rk_imex_5
Add a fifth order rk imex
2 parents 4d77602 + 12f38af commit bc1cb58

File tree

4 files changed

+205
-11
lines changed

4 files changed

+205
-11
lines changed

‎examples/acoustic_1d_imex/plot_dispersion.py‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ def findomega(stab_fh):
3030
sol = sol1
3131
else:
3232
print "Two roots with real part of same sign..."
33+
sol = sol0
3334
return sol
3435

3536
if __name__ == "__main__":
@@ -42,7 +43,7 @@ def findomega(stab_fh):
4243
swparams = {}
4344
swparams['collocation_class'] = collclass.CollGaussLegendre
4445
swparams['num_nodes'] = 3
45-
K = 4
46+
K = 5
4647
dirk_order = K
4748

4849
c_speed = 1.0
@@ -99,7 +100,7 @@ def findomega(stab_fh):
99100
# For testing, insert exact stability function exp(-dt*i*k*(Cs+Uadv)
100101
#stab_fh = la.expm(Cs+Uadv)
101102

102-
dirkts = dirk(Cs+Uadv, dirk_order)
103+
dirkts = dirk(Cs+Uadv, np.min([4,dirk_order]))
103104
stab_fh1 = dirkts.timestep(y1, 1.0)
104105
stab_fh2 = dirkts.timestep(y2, 1.0)
105106
stab_dirk = np.column_stack((stab_fh1, stab_fh2))

‎examples/acoustic_1d_imex/standard_integrators.py‎

Lines changed: 98 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
import math
33
import scipy.sparse.linalg as LA
44
import scipy.sparse as sp
5+
from decimal import *
56

67
#
78
# Runge-Kutta IMEX methods of order 1 to 3
@@ -13,7 +14,7 @@ def __init__(self, M_fast, M_slow, order):
1314
assert np.shape(M_slow)[0]==np.shape(M_slow)[1], "A_slow must be square"
1415
assert np.shape(M_fast)[0]==np.shape(M_slow)[0], "A_fast and A_slow must be of the same size"
1516

16-
assert order in [1,2,3,4], "Order must be 1, 2, 3 or 4"
17+
assert order in [1,2,3,4,5], "Order must be between 1 and 5"
1718
self.order = order
1819

1920
if self.order==1:
@@ -58,7 +59,102 @@ def __init__(self, M_fast, M_slow, order):
5859
self.b = np.array([82889./524892.,0,15625./83664.,69875./102672.,-2260./8211,1./4])
5960
self.b_hat = np.array([4586570599./29645900160.,0,178811875./945068544.,814220225./1159782912.,-3700637./11593932.,61727./225920.])
6061
self.nstages = 6
61-
62+
63+
elif self.order==5:
64+
65+
# from Kennedy and Carpenter
66+
# copied from http://www.mcs.anl.gov/petsc/petsc-3.2/src/ts/impls/arkimex/arkimex.c
67+
self.A_hat = np.zeros((8,8))
68+
getcontext().prec = 56
69+
self.A_hat[1,0] = Decimal(41.0)/Decimal(100.0)
70+
self.A_hat[2,0] = Decimal(367902744464.)/Decimal(2072280473677.)
71+
self.A_hat[2,1] = Decimal(677623207551.)/Decimal(8224143866563.)
72+
self.A_hat[3,0] = Decimal(1268023523408.)/Decimal(10340822734521.)
73+
self.A_hat[3,1] = 0.0
74+
self.A_hat[3,2] = Decimal(1029933939417.)/Decimal(13636558850479.)
75+
self.A_hat[4,0] = Decimal(14463281900351.)/Decimal(6315353703477.)
76+
self.A_hat[4,1] = 0.0
77+
self.A_hat[4,2] = Decimal(66114435211212.)/Decimal(5879490589093.)
78+
self.A_hat[4,3] = Decimal(-54053170152839.)/Decimal(4284798021562.)
79+
self.A_hat[5,0] = Decimal(14090043504691.)/Decimal(34967701212078.)
80+
self.A_hat[5,1] = 0.0
81+
self.A_hat[5,2] = Decimal(15191511035443.)/Decimal(11219624916014.)
82+
self.A_hat[5,3] = Decimal(-18461159152457.)/Decimal(12425892160975.)
83+
self.A_hat[5,4] = Decimal(-281667163811.)/Decimal(9011619295870.)
84+
self.A_hat[6,0] = Decimal(19230459214898.)/Decimal(13134317526959.)
85+
self.A_hat[6,1] = 0.0
86+
self.A_hat[6,2] = Decimal(21275331358303.)/Decimal(2942455364971.)
87+
self.A_hat[6,3] = Decimal(-38145345988419.)/Decimal(4862620318723.)
88+
self.A_hat[6,4] = Decimal(-1.0)/Decimal(8.0)
89+
self.A_hat[6,5] = Decimal(-1.0)/Decimal(8.0)
90+
self.A_hat[7,0] = Decimal(-19977161125411.)/Decimal(11928030595625.)
91+
self.A_hat[7,1] = 0.0
92+
self.A_hat[7,2] = Decimal(-40795976796054.)/Decimal(6384907823539.)
93+
self.A_hat[7,3] = Decimal(177454434618887.)/Decimal(12078138498510.)
94+
self.A_hat[7,4] = Decimal(782672205425.)/Decimal(8267701900261.)
95+
self.A_hat[7,5] = Decimal(-69563011059811.)/Decimal(9646580694205.)
96+
self.A_hat[7,6] = Decimal(7356628210526.)/Decimal(4942186776405.)
97+
98+
self.b_hat = np.zeros(8)
99+
self.b_hat[0] = Decimal(-872700587467.)/Decimal(9133579230613.)
100+
self.b_hat[1] = 0.0
101+
self.b_hat[2] = 0.0
102+
self.b_hat[3] = Decimal(22348218063261.)/Decimal(9555858737531.)
103+
self.b_hat[4] = Decimal(-1143369518992.)/Decimal(8141816002931.)
104+
self.b_hat[5] = Decimal(-39379526789629.)/Decimal(19018526304540.)
105+
self.b_hat[6] = Decimal(32727382324388.)/Decimal(42900044865799.)
106+
self.b_hat[7] = Decimal(41.0)/Decimal(200.0)
107+
108+
self.A = np.zeros((8,8))
109+
self.A[1,0] = Decimal(41.)/Decimal(200.)
110+
self.A[1,1] = Decimal(41.)/Decimal(200.)
111+
self.A[2,0] = Decimal(41.)/Decimal(400.)
112+
self.A[2,1] = Decimal(-567603406766.)/Decimal(11931857230679.)
113+
self.A[2,2] = Decimal(41.)/Decimal(200.)
114+
self.A[3,0] = Decimal(683785636431.)/Decimal(9252920307686.)
115+
self.A[3,1] = 0.0
116+
self.A[3,2] = Decimal(-110385047103.)/Decimal(1367015193373.)
117+
self.A[3,3] = Decimal(41.)/Decimal(200.)
118+
self.A[4,0] = Decimal(3016520224154.)/Decimal(10081342136671.)
119+
self.A[4,1] = 0.0
120+
self.A[4,2] = Decimal(30586259806659.)/Decimal(12414158314087.)
121+
self.A[4,3] = Decimal(-22760509404356.)/Decimal(11113319521817.)
122+
self.A[4,4] = Decimal(41.)/Decimal(200.)
123+
self.A[5,0] = Decimal(218866479029.)/Decimal(1489978393911.)
124+
self.A[5,1] = 0.0
125+
self.A[5,2] = Decimal(638256894668.)/Decimal(5436446318841.)
126+
self.A[5,3] = Decimal(-1179710474555.)/Decimal(5321154724896.)
127+
self.A[5,4] = Decimal(-60928119172.)/Decimal(8023461067671.)
128+
self.A[5,5] = Decimal(41.)/Decimal(200.)
129+
self.A[6,0] = Decimal(1020004230633.)/Decimal(5715676835656.)
130+
self.A[6,1] = 0.0
131+
self.A[6,2] = Decimal(25762820946817.)/Decimal(25263940353407.)
132+
self.A[6,3] = Decimal(-2161375909145.)/Decimal(9755907335909.)
133+
self.A[6,4] = Decimal(-211217309593.)/Decimal(5846859502534.)
134+
self.A[6,5] = Decimal(-4269925059573.)/Decimal(7827059040749.)
135+
self.A[6,6] = Decimal(41.)/Decimal(200.)
136+
self.A[7,0] = Decimal(-872700587467.)/Decimal(9133579230613.)
137+
self.A[7,1] = 0.0
138+
self.A[7,2] = 0.0
139+
self.A[7,3] = Decimal(22348218063261.)/Decimal(9555858737531.)
140+
self.A[7,4] = Decimal(-1143369518992.)/Decimal(8141816002931.)
141+
self.A[7,5] = Decimal(-39379526789629.)/Decimal(19018526304540.)
142+
self.A[7,6] = Decimal(32727382324388.)/Decimal(42900044865799.)
143+
self.A[7,7] = Decimal(41.)/Decimal(200.)
144+
145+
self.b = np.zeros(8)
146+
147+
self.b[0] = Decimal(-975461918565.)/Decimal(9796059967033.)
148+
self.b[1] = 0.0
149+
self.b[2] = 0.0
150+
self.b[3] = Decimal(78070527104295.)/Decimal(32432590147079.)
151+
self.b[4] = Decimal(-548382580838.)/Decimal(3424219808633.)
152+
self.b[5] = Decimal(-33438840321285.)/Decimal(15594753105479.)
153+
self.b[6] = Decimal(3629800801594.)/Decimal(4656183773603.)
154+
self.b[7] = Decimal(4035322873751.)/Decimal(18575991585200.)
155+
156+
self.nstages = 8
157+
62158
self.M_fast = sp.csc_matrix(M_fast)
63159
self.M_slow = sp.csc_matrix(M_slow)
64160
self.ndof = np.shape(M_fast)[0]

‎examples/boussinesq_2d_imex/rungmrescounter.py‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -40,14 +40,14 @@
4040
swparams['do_LU'] = False
4141

4242
sparams = {}
43-
sparams['maxiter'] = 3
43+
sparams['maxiter'] = 5
4444

45-
dirk_order = 3
45+
dirk_order = 5
4646

4747
# setup parameters "in time"
4848
t0 = 0
4949
Tend = 3000
50-
Nsteps = 100
50+
Nsteps = 500
5151
dt = Tend/float(Nsteps)
5252

5353
# This comes as read-in for the problem class
@@ -63,7 +63,7 @@
6363
pparams['gmres_maxiter'] = [500]
6464
pparams['gmres_restart'] = [10]
6565
pparams['gmres_tol_limit'] = [1e-5]
66-
pparams['gmres_tol_factor'] = [0.05]
66+
pparams['gmres_tol_factor'] = [0.1]
6767

6868
# This comes as read-in for the transfer operations
6969
tparams = {}
@@ -96,7 +96,7 @@
9696
print("CFL number of acoustics (horizontal): %4.2f" % cfl_acoustic_hor)
9797
print("CFL number of acoustics (vertical): %4.2f" % cfl_acoustic_ver)
9898

99-
dirkp = dirk(P, dirk_order)
99+
dirkp = dirk(P, np.min([4,dirk_order]))
100100
u0 = uinit.values.flatten()
101101
udirk = np.copy(u0)
102102
print "Running DIRK ...."
@@ -116,7 +116,7 @@
116116

117117
# For reference solution, increase GMRES tolerance
118118
P.gmres_tol_limit = 1e-10
119-
rkimexref = rk_imex(P, 4)
119+
rkimexref = rk_imex(P, 5)
120120
uref = np.copy(u0)
121121
dt_ref = dt/10.0
122122
print "Running RK-IMEX reference...."

‎examples/boussinesq_2d_imex/standard_integrators.py‎

Lines changed: 98 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
import scipy.sparse.linalg as LA
44
import scipy.sparse as sp
55
from ProblemClass import Callback, logging, boussinesq_2d_imex
6+
from decimal import *
67

78
#
89
# Runge-Kutta IMEX methods of order 1 to 3
@@ -11,7 +12,7 @@ class rk_imex:
1112

1213
def __init__(self, problem, order):
1314

14-
assert order in [1,2,3,4], "Order must be 1, 2, 3 or 4"
15+
assert order in [1,2,3,4,5], "Order must be between 1 and 5"
1516
self.order = order
1617

1718
if self.order==1:
@@ -57,6 +58,102 @@ def __init__(self, problem, order):
5758
self.b_hat = np.array([4586570599./29645900160.,0,178811875./945068544.,814220225./1159782912.,-3700637./11593932.,61727./225920.])
5859
self.nstages = 6
5960

61+
elif self.order==5:
62+
63+
# from Kennedy and Carpenter
64+
# copied from http://www.mcs.anl.gov/petsc/petsc-3.2/src/ts/impls/arkimex/arkimex.c
65+
self.A_hat = np.zeros((8,8))
66+
getcontext().prec = 56
67+
self.A_hat[1,0] = Decimal(41.0)/Decimal(100.0)
68+
self.A_hat[2,0] = Decimal(367902744464.)/Decimal(2072280473677.)
69+
self.A_hat[2,1] = Decimal(677623207551.)/Decimal(8224143866563.)
70+
self.A_hat[3,0] = Decimal(1268023523408.)/Decimal(10340822734521.)
71+
self.A_hat[3,1] = 0.0
72+
self.A_hat[3,2] = Decimal(1029933939417.)/Decimal(13636558850479.)
73+
self.A_hat[4,0] = Decimal(14463281900351.)/Decimal(6315353703477.)
74+
self.A_hat[4,1] = 0.0
75+
self.A_hat[4,2] = Decimal(66114435211212.)/Decimal(5879490589093.)
76+
self.A_hat[4,3] = Decimal(-54053170152839.)/Decimal(4284798021562.)
77+
self.A_hat[5,0] = Decimal(14090043504691.)/Decimal(34967701212078.)
78+
self.A_hat[5,1] = 0.0
79+
self.A_hat[5,2] = Decimal(15191511035443.)/Decimal(11219624916014.)
80+
self.A_hat[5,3] = Decimal(-18461159152457.)/Decimal(12425892160975.)
81+
self.A_hat[5,4] = Decimal(-281667163811.)/Decimal(9011619295870.)
82+
self.A_hat[6,0] = Decimal(19230459214898.)/Decimal(13134317526959.)
83+
self.A_hat[6,1] = 0.0
84+
self.A_hat[6,2] = Decimal(21275331358303.)/Decimal(2942455364971.)
85+
self.A_hat[6,3] = Decimal(-38145345988419.)/Decimal(4862620318723.)
86+
self.A_hat[6,4] = Decimal(-1.0)/Decimal(8.0)
87+
self.A_hat[6,5] = Decimal(-1.0)/Decimal(8.0)
88+
self.A_hat[7,0] = Decimal(-19977161125411.)/Decimal(11928030595625.)
89+
self.A_hat[7,1] = 0.0
90+
self.A_hat[7,2] = Decimal(-40795976796054.)/Decimal(6384907823539.)
91+
self.A_hat[7,3] = Decimal(177454434618887.)/Decimal(12078138498510.)
92+
self.A_hat[7,4] = Decimal(782672205425.)/Decimal(8267701900261.)
93+
self.A_hat[7,5] = Decimal(-69563011059811.)/Decimal(9646580694205.)
94+
self.A_hat[7,6] = Decimal(7356628210526.)/Decimal(4942186776405.)
95+
96+
self.b_hat = np.zeros(8)
97+
self.b_hat[0] = Decimal(-872700587467.)/Decimal(9133579230613.)
98+
self.b_hat[1] = 0.0
99+
self.b_hat[2] = 0.0
100+
self.b_hat[3] = Decimal(22348218063261.)/Decimal(9555858737531.)
101+
self.b_hat[4] = Decimal(-1143369518992.)/Decimal(8141816002931.)
102+
self.b_hat[5] = Decimal(-39379526789629.)/Decimal(19018526304540.)
103+
self.b_hat[6] = Decimal(32727382324388.)/Decimal(42900044865799.)
104+
self.b_hat[7] = Decimal(41.0)/Decimal(200.0)
105+
106+
self.A = np.zeros((8,8))
107+
self.A[1,0] = Decimal(41.)/Decimal(200.)
108+
self.A[1,1] = Decimal(41.)/Decimal(200.)
109+
self.A[2,0] = Decimal(41.)/Decimal(400.)
110+
self.A[2,1] = Decimal(-567603406766.)/Decimal(11931857230679.)
111+
self.A[2,2] = Decimal(41.)/Decimal(200.)
112+
self.A[3,0] = Decimal(683785636431.)/Decimal(9252920307686.)
113+
self.A[3,1] = 0.0
114+
self.A[3,2] = Decimal(-110385047103.)/Decimal(1367015193373.)
115+
self.A[3,3] = Decimal(41.)/Decimal(200.)
116+
self.A[4,0] = Decimal(3016520224154.)/Decimal(10081342136671.)
117+
self.A[4,1] = 0.0
118+
self.A[4,2] = Decimal(30586259806659.)/Decimal(12414158314087.)
119+
self.A[4,3] = Decimal(-22760509404356.)/Decimal(11113319521817.)
120+
self.A[4,4] = Decimal(41.)/Decimal(200.)
121+
self.A[5,0] = Decimal(218866479029.)/Decimal(1489978393911.)
122+
self.A[5,1] = 0.0
123+
self.A[5,2] = Decimal(638256894668.)/Decimal(5436446318841.)
124+
self.A[5,3] = Decimal(-1179710474555.)/Decimal(5321154724896.)
125+
self.A[5,4] = Decimal(-60928119172.)/Decimal(8023461067671.)
126+
self.A[5,5] = Decimal(41.)/Decimal(200.)
127+
self.A[6,0] = Decimal(1020004230633.)/Decimal(5715676835656.)
128+
self.A[6,1] = 0.0
129+
self.A[6,2] = Decimal(25762820946817.)/Decimal(25263940353407.)
130+
self.A[6,3] = Decimal(-2161375909145.)/Decimal(9755907335909.)
131+
self.A[6,4] = Decimal(-211217309593.)/Decimal(5846859502534.)
132+
self.A[6,5] = Decimal(-4269925059573.)/Decimal(7827059040749.)
133+
self.A[6,6] = Decimal(41.)/Decimal(200.)
134+
self.A[7,0] = Decimal(-872700587467.)/Decimal(9133579230613.)
135+
self.A[7,1] = 0.0
136+
self.A[7,2] = 0.0
137+
self.A[7,3] = Decimal(22348218063261.)/Decimal(9555858737531.)
138+
self.A[7,4] = Decimal(-1143369518992.)/Decimal(8141816002931.)
139+
self.A[7,5] = Decimal(-39379526789629.)/Decimal(19018526304540.)
140+
self.A[7,6] = Decimal(32727382324388.)/Decimal(42900044865799.)
141+
self.A[7,7] = Decimal(41.)/Decimal(200.)
142+
143+
self.b = np.zeros(8)
144+
145+
self.b[0] = Decimal(-975461918565.)/Decimal(9796059967033.)
146+
self.b[1] = 0.0
147+
self.b[2] = 0.0
148+
self.b[3] = Decimal(78070527104295.)/Decimal(32432590147079.)
149+
self.b[4] = Decimal(-548382580838.)/Decimal(3424219808633.)
150+
self.b[5] = Decimal(-33438840321285.)/Decimal(15594753105479.)
151+
self.b[6] = Decimal(3629800801594.)/Decimal(4656183773603.)
152+
self.b[7] = Decimal(4035322873751.)/Decimal(18575991585200.)
153+
154+
self.nstages = 8
155+
156+
60157
self.problem = problem
61158
self.ndof = np.shape(problem.M)[0]
62159
self.logger = logging()

0 commit comments

Comments
 (0)