-
Notifications
You must be signed in to change notification settings - Fork 0
/
gark_methods.py
59 lines (51 loc) · 2.4 KB
/
gark_methods.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
import numpy as np
from additive_rk import ark_solve
from butcher_tableau import Tableau
def gark_convert(A, b):
""" Convert the provided GARK structure to an ARK structure """
# This assumes M = N, and J = I
tableaus = []
Ni = sum([x.size for x in b])
pre = 0
for i in range(len(b)):
si = b[i].size
ai = np.zeros((Ni, Ni))
ai[:,pre:pre+si] = np.concatenate([A[q][i] for q in range(len(b))])
bi = np.zeros(Ni)
bi[pre:pre+si] = b[i]
ci = ai.sum(axis=1)
tableaus.append(Tableau(ci, ai, bi))
pre += si
return tableaus
def gark_solve(f, dt, y0, t0, tf, A, b, bc=None, solver_parameters={}, fname=None, save_steps = 0, jacobian=None):
""" This function uses an GARK method to solve a differential equation
This is done by converting to an ARK structure
-----
Inputs:
functions - the operators of the differential equation, each with inputs (t, y)
these may also be finite element Forms as provided by firedrake
dt - the amount time will increase by
y0 - the current value of y
if using the finite element version, this should be of type Function
t0 - the current value of t
if using the finite element version, this should be of type Constant
tf - the time to solve until
A - the list of lists containing the component arrays a{J(q), i}
b - the list containing the b{i} arrays
bc - optional. Only used if using the finite element version. The boundary condition(s) to apply.
solver_parameters - optional. Only used for the finite element version. Any solver parameters to use (see firedrake documentation for details)
fname - optional. If provided, will save intermediate results to this file.
- if using the finite element version of the code, this is a HDF5 file
otherwise it is a csv file.
save_steps - the number of intermediate steps to save if fname is provided
if it is not provided, the default is to save every step
(or after every dt if embedded methods are being used).
Return
-----
the approximate value of y at tf
"""
if len(f) > len(b):
print("ERROR: not enough tableau provided")
return
tableaus = gark_convert(A, b)
return ark_solve(f, dt, y0, t0, tf, tableaus, bc=bc, solver_parameters=solver_parameters, fname=fname, save_steps = save_steps, jacobian=jacobian)