/[escript]/trunk/ripley/test/python/cook1.py
ViewVC logotype

Contents of /trunk/ripley/test/python/cook1.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3479 - (show annotations)
Wed Mar 23 04:03:22 2011 UTC (8 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 6382 byte(s)


1 import numpy as np
2 import matplotlib.pyplot as plt
3 import random
4
5 from esys.escript import *
6 from esys.escript.linearPDEs import LinearPDE
7 from esys.escript.pdetools import Locator
8 from esys.finley import Rectangle,Brick
9
10 import scipy.optimize as so
11
12 NE=200
13 NE=100
14 H=1.
15 padding=0.0
16 numContours=25
17 numPlotsGrid=25
18
19 class PoissonTest(object):
20 def __init__(self, domain, a, b, xc):
21 self.domain=domain
22 self.a=a
23 self.b=b
24 self.xc=xc
25 self.pde=LinearPDE(domain)
26 self.pde.setSymmetryOn()
27 x=self.domain.getX()
28 DIM=domain.getDim()
29 self.pde.setValue(A=kronecker(DIM), q=whereZero(x[DIM-1]-inf(x[DIM-1])))
30 self.u_ref=self.getU(a, b, xc)
31 self.omega=1.
32
33 def getU(self, a, b, xc):
34 x=Function(self.domain).getX()
35 self.pde.setValue(Y=exp(a-(b*length(x-xc))**2))
36 u=self.pde.getSolution()
37 print "PoissonTest.getU: a,b,xc, u =",a,b,xc,inf(u),sup(u)
38 return u
39
40 def getCostFunction(self, a, b, xc):
41 u=interpolate(self.getU(a, b, xc), FunctionOnBoundary(self.domain))
42 C=0.5*integrate( self.omega * (u-self.u_ref)**2)
43 return C
44
45 class PoissonTestTopData(PoissonTest):
46 def __init__(self, domain, a, b, xc):
47 PoissonTest.__init__(self, domain, a, b, xc)
48 x=FunctionOnBoundary(self.domain).getX()
49 DIM=self.domain.getDim()
50 self.omega=whereZero(x[DIM-1]-sup(x[DIM-1]))
51 class PoissonTestDiscreteData(PoissonTest):
52 def __init__(self, N, sigma, domain, a, b, xc):
53 PoissonTest.__init__(self, domain, a, b, xc)
54 self.N=N
55 self.sigma=sigma
56 x=domain.getX()
57 x_min=inf(x[0])
58 x_max=sup(x[0])
59 z_max=sup(x[self.domain.getDim()-1])
60 dx=(x_max-x_min)/N
61 self.loc=Locator(self.domain, [ (dx/2+dx*i, z_max) for i in xrange(N) ])
62 uu=self.loc(self.u_ref)
63 s=self.sigma*(max(uu)+min(uu))/2.
64 print self.loc.getX()
65 self.data = numpy.array([ random.gauss(d, s) for d in self.loc(self.u_ref) ])
66 print self.data
67 def getCostFunction(self, a, b, xc):
68 u=numpy.array(self.loc(self.getU(a, b, xc)))
69 C=(0.5/self.N)*numpy.sum((u-self.data)**2)
70 return C
71
72
73
74 def plotContour(p, filename, title):
75 x0=np.linspace(-padding, 1+padding, numPlotsGrid)
76 X0, Y0 = np.meshgrid(x0, H*x0)
77
78 plt.figure()
79 C=X0.copy()
80 for i in range(C.shape[0]):
81 for j in range(C.shape[1]):
82 C[i,j]=p.getCostFunction(p.a,p.b, (X0[i,j], Y0[i,j]))
83
84 C_max=np.amax(C)
85 C_min=np.amin(C)
86 levels = np.linspace(C_min, C_max, numContours)
87 cset1 = plt.contourf(X0, Y0, C, levels, cmap=plt.cm.get_cmap('jet', len(levels)-1))
88 cset2 = plt.contour(X0, Y0, C, cset1.levels, colors = 'k', hold='on')
89 for c in cset2.collections: c.set_linestyle('solid')
90 plt.title('Cost function (%s)'%title)
91 plt.xlabel("xc_0")
92 plt.ylabel("xc_1")
93 plt.colorbar(cset1)
94 plt.savefig(filename)
95
96
97
98 plotContour(PoissonTestDiscreteData(5, 0.10, Rectangle(NE, int(NE*H),l1=H), 0,15,(0.5,0.5)), "cost_passion_discrete_w_cc.png","b=10, xc=(0.5,0.2)")
99
100 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.2)), "cost_passion_top_w_cb.png","b=10, xc=(0.5,0.2)")
101 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.5)), "cost_passion_top_w_cc.png","b=10, xc=(0.5,0.5)")
102 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.8)), "cost_passion_top_w_ct.png","b=10, xc=(0.5,0.8)")
103 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.2)), "cost_passion_top_w_lb.png","b=10, xc=(0.2,0.2)")
104 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.5)), "cost_passion_top_w_lc.png","b=10, xc=(0.2,0.5)")
105 # plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.8)), "cost_passion_top_w_lt.png","b=10, xc=(0.2,0.8)")
106 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.2)), "cost_passion_all_w_cb.png","b=10, xc=(0.5,0.2)")
107 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.5)), "cost_passion_all_w_cc.png","b=10, xc=(0.5,0.5)")
108 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.5,0.8)), "cost_passion_all_w_ct.png","b=10, xc=(0.5,0.8)")
109 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.2)), "cost_passion_all_w_lb.png","b=10, xc=(0.2,0.2)")
110 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.5)), "cost_passion_all_w_lc.png","b=10, xc=(0.2,0.5)")
111 # plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,10,(0.2,0.8)), "cost_passion_all_w_lt.png","b=10, xc=(0.2,0.8)")
112
113 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.2)), "cost_passion_top_n_cb.png","b=100, xc=(0.5,0.2)")
114 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.5)), "cost_passion_top_n_cc.png","b=100, xc=(0.5,0.5)")
115 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.8)), "cost_passion_top_n_ct.png","b=100, xc=(0.5,0.8)")
116 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.2)), "cost_passion_top_n_lb.png","b=100, xc=(0.2,0.2)")
117 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.5)), "cost_passion_top_n_lc.png","b=100, xc=(0.2,0.5)")
118 #plotContour(PoissonTestTopData(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.8)), "cost_passion_top_n_lt.png","b=100, xc=(0.2,0.8)")
119 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.2)), "cost_passion_all_n_cb.png","b=100, xc=(0.5,0.2)")
120 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.5)), "cost_passion_all_n_cc.png","b=100, xc=(0.5,0.5)")
121 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.5,0.8)), "cost_passion_all_n_ct.png","b=100, xc=(0.5,0.8)")
122 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.2)), "cost_passion_all_n_lb.png","b=100, xc=(0.2,0.2)")
123 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.5)), "cost_passion_all_n_lc.png","b=100, xc=(0.2,0.5)")
124 #plotContour(PoissonTest(Rectangle(NE, int(NE*H),l1=H), 0,100,(0.2,0.8)), "cost_passion_all_n_lt.png","b=100, xc=(0.2,0.8)")
125
126
127
128
129

  ViewVC Help
Powered by ViewVC 1.1.26