-
Notifications
You must be signed in to change notification settings - Fork 2
/
python_kernel_test.py
175 lines (139 loc) · 6.37 KB
/
python_kernel_test.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
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
# -*- coding: utf-8 -*-
#
# Tests/usage examples for a custom Python-based kernel.
from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt
from pydgq.solver.types import DTYPE
from pydgq.solver.galerkin import init
from pydgq.solver.kernel_interface import PythonKernel
import pydgq.solver.odesolve
from pydgq.utils.discontify import discontify # for plotting dG results
#####################
# config for testing
#####################
q = 2 # degree of basis for dG and cG
# How many visualization (interpolation) points to use within each timestep for Galerkin methods.
#
# Note that the dG solution has the best accuracy at the endpoint of the timestep;
# to compare apples-to-apples with classical integrators, this should be set to 1.
#
# Larger values (e.g. 11) are useful for visualizing the behavior of the dG solution inside
# the timestep (something the classical integrators do not model at all).
#
nt_vis_galerkin = 1
nt = 100 # number of timesteps
dt = 0.1 # timestep size
save_from = 0 # see pydgq.solver.odesolve.ivp()
#####################
# custom kernel
#####################
# A simple cosine kernel with phase-shifted components.
#
# See cython_kernel.pyx for a Cython-accelerated version of this.
#
# The custom kernel only needs to override callback(); even __init__ is not strictly needed,
# unless adding some custom parameters (like here).
#
class MyKernel(PythonKernel):
def __init__(self, n, omega): # omega : rad/s
# super
PythonKernel.__init__(self, n)
# custom init
self.omega = omega
def callback(self, t):
for j in range(self.n):
phi0_j = (float(j+1) / self.n) * 2. * np.pi
self.out[j] = np.cos(phi0_j + self.omega*t)
# known analytical solution, for testing the integrators
#
def reference_solution(self, tt):
tt = np.atleast_1d(tt)
ww = np.empty( (tt.shape[0],self.n), dtype=DTYPE, order="C" )
sol = lambda t,phi0 : 1./self.omega * np.sin(phi0 + self.omega*t)
for j in range(self.n):
phi0_j = (float(j+1) / self.n) * 2. * np.pi
ww[:,j] = sol(tt,phi0_j) - sol(0.,phi0_j) # shift to account for the initial condition (all solution components start at zero)
return ww
#####################
# main program
#####################
# rel_tol : how close the numerical solution must be to the exact analytical one, in l-infinity norm (max abs)
def test(integrator, nt_vis, rel_tol=1e-2, vis=False):
n_saved_timesteps = pydgq.solver.odesolve.n_saved_timesteps( nt, save_from )
result_len = pydgq.solver.odesolve.result_len( nt, save_from, interp=nt_vis )
startj,endj = pydgq.solver.odesolve.timestep_boundaries( nt, save_from, interp=nt_vis )
n = 3 # number of DOFs in the 1st-order system
w0 = np.zeros( (n,), dtype=DTYPE, order="C" ) # a trivial IC
# instantiate kernel
rhs = MyKernel(n, omega=0.1 * (2. * np.pi))
# create output arrays
ww = None #np.empty( (result_len,n), dtype=DTYPE, order="C" ) # result array for w; if None, will be created by ivp()
ff = np.empty( (result_len,n), dtype=DTYPE, order="C" ) # optional, result array for w', could be None
fail = np.empty( (n_saved_timesteps,), dtype=np.intc, order="C" ) # optional, fail flag for each timestep, could be None
# solve problem
ww,tt = pydgq.solver.odesolve.ivp( integrator=integrator, allow_denormals=False,
w0=w0, dt=dt, nt=nt,
save_from=save_from, interp=nt_vis,
rhs=rhs,
ww=ww, ff=ff, fail=fail,
maxit=100 )
# check result
ww_ref = rhs.reference_solution(tt)
relerr_linfty = np.linalg.norm(ww - ww_ref, ord=np.inf) / np.linalg.norm(ww_ref, ord=np.inf)
if (relerr_linfty < rel_tol).all():
passed = 1
if not vis:
print("PASS, %03s, tol=% 7g, relerr=%g" % (integrator, rel_tol, relerr_linfty))
else:
passed = 0
if not vis:
print("FAIL, %03s, tol=% 7g, relerr=%g" % (integrator, rel_tol, relerr_linfty))
# visualize if requested
if vis:
plt.figure(1)
plt.clf()
# show the discontinuities at timestep boundaries if using dG (and actually have something to draw within each timestep)
if integrator == "dG" and nt_vis > 1:
tt = discontify( tt, endj - 1, fill="nan" )
for j in range(n):
# we need the copy() to get memory-contiguous data for discontify() to process
wtmp = discontify( ww_ref[:,j].copy(), endj - 1, fill="nan" )
plt.plot( tt, wtmp, 'r--' ) # known exact solution
wtmp = discontify( ww[:,j].copy(), endj - 1, fill="nan" )
plt.plot( tt, wtmp, 'k-' ) # numerical solution
else:
plt.plot( tt, ww_ref, 'r--' ) # known exact solution
plt.plot( tt, ww, 'k-' ) # numerical solution
plt.grid(b=True, which="both")
plt.axis("tight")
plt.xlabel(r"$t$")
plt.ylabel(r"$w(t)$")
return passed
if __name__ == '__main__':
print("** Testing integrators with custom Python kernel **")
# "SE" is not applicable, since we are testing a 1st-order problem
stuff_to_test = ( ("IMR", 1e-3, False),
("BE", 1e-1, False),
("RK4", 1e-8, False),
("RK3", 1e-6, False),
("RK2", 1e-3, False),
("FE", 1e-1, False),
("dG", 1e-12, True ),
("cG", 1e-4, True )
)
n_passed = 0
for integrator,rel_tol,is_galerkin in stuff_to_test:
if is_galerkin:
nt_vis = nt_vis_galerkin # visualization points per timestep in Galerkin methods
init(q=q, method=integrator, nt_vis=nt_vis, rule=None)
else:
nt_vis = 1 # other methods compute only the end value for each timestep
n_passed += test(integrator, nt_vis, rel_tol)
print("** %d/%d tests passed **" % (n_passed, len(stuff_to_test)))
# # DEBUG: draw something
# #
# nt_vis = nt_vis_galerkin
# init(q=q, method="dG", nt_vis=nt_vis, rule=None)
# test(integrator="dG", nt_vis=nt_vis, vis=True)
# plt.show()