made a get_error function that is generic and functional. I'm still getting unsavory numbers for the cubic, 2-D case on my single test case. I will have to try out more refined grids, albeit the quadratic is consistently improving results for my case.

This commit is contained in:
Stephen McQuay 2010-05-05 23:03:12 -06:00
parent 3a1c13bcac
commit df35cb174b
4 changed files with 59 additions and 12 deletions

View File

@ -1,4 +1,4 @@
#!/usr/bin/python #!/usr/bin/env python
import sys import sys
from optparse import OptionParser from optparse import OptionParser

View File

@ -37,7 +37,7 @@ if __name__ == '__main__':
sys.exit(1) sys.exit(1)
print "total points", total_points smblog.debug("total points: %d" % total_points)
@ -67,3 +67,4 @@ if __name__ == '__main__':
good.append(d[True]) good.append(d[True])
draw_gb(bad = bad, good = good) draw_gb(bad = bad, good = good)

View File

@ -1,6 +1,5 @@
#!/usr/bin/env python #!/usr/bin/env python
import math
from baker import run_baker from baker import run_baker
from baker.tools import improved_answer from baker.tools import improved_answer
@ -9,10 +8,12 @@ from baker.tools import smblog
from grid.DD import grid from grid.DD import grid
from grid.simplex import contains from grid.simplex import contains
import numpy as np
def exact_func(X): def exact_func(X):
x = X[0] x = X[0]
y = X[0] y = X[0]
return 1 - math.sin((x-0.5)**2 + (y-0.5)**2) return 1 - np.sin(np.power((x-0.5), 2) + np.power((y-0.5), 2))
if __name__ == '__main__': if __name__ == '__main__':
points = [ points = [
@ -33,8 +34,6 @@ if __name__ == '__main__':
g.construct_connectivity() g.construct_connectivity()
R = g.create_mesh(range(3)) R = g.create_mesh(range(3))
print contains(X, R.points)
try: try:
extra = int(sys.argv[1]) extra = int(sys.argv[1])
except: except:

View File

@ -88,7 +88,7 @@ def qlinear(X, R):
""" """
phis = get_phis(X, R.points) phis = get_phis(X, R.points)
qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)]) qlin = np.sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)])
return phis, qlin return phis, qlin
def qlinear_3D(X, R): def qlinear_3D(X, R):
@ -185,6 +185,54 @@ def get_error_cubic(phi, R, S):
return error_term return error_term
def get_error_sauron(phi, R, S, order = 2):
smblog.debug("len(phi): %d"% len(phi))
B = [] # baker eq 9
w = [] # baker eq 11
p = pattern(order, len(phi), offset = -1)
smblog.debug("pattern: %s" % p)
for (s,q) in zip(S.points, S.q):
cur_phi, cur_qlin = qlinear(s, R)
l = []
for i in p:
cur_sum = cur_phi[i[0]]
for j in i[1:]:
cur_sum *= cur_phi[j]
l.append(cur_sum)
B.append(l)
w.append(q - cur_qlin)
smblog.debug("B: %s" % B)
smblog.debug("w: %s" % w)
B = np.array(B)
w = np.array(w)
A = np.dot(B.T, B)
b = np.dot(B.T, w)
# baker solve eq 10
try:
abc = np.linalg.solve(A,b)
except np.linalg.LinAlgError as e:
smblog.error("linear calculation went bad, resorting to np.linalg.pinv: %s" % e)
abc = np.dot(np.linalg.pinv(A), b)
smblog.debug(len(abc) == len(p))
error_term = 0.0
for (a, i) in zip(abc, p):
cur_sum = a
for j in i:
cur_sum *= phi[j]
error_term += cur_sum
smblog.debug("error_term smb: %s" % error_term)
return error_term, abc
def run_baker(X, R, S, order=2): def run_baker(X, R, S, order=2):
""" """
@ -196,6 +244,7 @@ def run_baker(X, R, S, order=2):
R = Simplex R = Simplex
S = extra points S = extra points
""" """
smblog.debug("order = %d" % order)
answer = { answer = {
'qlin': None, 'qlin': None,
@ -208,12 +257,10 @@ def run_baker(X, R, S, order=2):
if order == 1: if order == 1:
answer['qlin'] = qlin answer['qlin'] = qlin
return answer return answer
elif order == 2: elif order in (2,3):
error_term = get_error_quadratic(phi, R, S) error_term, abc = get_error_sauron(phi, R, S, order)
elif order == 3:
error_term = get_error_cubic(phi, R, S)
else: else:
raise smberror('unacceptable order for baker method') raise smberror('unsupported order for baker method')
q_final = qlin + error_term q_final = qlin + error_term