test/python/SimpleSolve.py

# $Id$

import sys
import os
import unittest

from esys.escript import *
from esys.linearPDEs import *
from esys import finley

print "\nSimpleSolve.py"
print "--------------"

alpha=0.01

# generate mesh

# print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
# mydomain=finley.Rectangle(1,1)

# print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
mydomain=finley.Rectangle(4,4,1)

print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
# mydomain=finley.Rectangle(151,151,1)

print "\nSetup domain and functions"
print "--------------------------"

print "e=Function(mydomain):"
e=Function(mydomain)

print "n=ContinuousFunction(mydomain):"
n=ContinuousFunction(mydomain)

# get handles to nodes and elements 1

print "\nGet handles to nodes and elements(1)=>"
print "--------------------------------------"

print "u_ex=Scalar(1,n,True):"
u_ex=Scalar(1,n,True)

print "x=e.getX():"
x=e.getX()

print "norm_u_ex=u_ex.Lsup():"
norm_u_ex=u_ex.Lsup()

print "\nGenerate a test solution (1)"
print "----------------------------"

print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
mypde=LinearPDE(mydomain)
mypde.setDebugOn()
mypde.setValue(A=[[1.,0.8],[0.4,1.]],D=alpha,Y=alpha)

print "mypde.checkSymmetry()"
print mypde.checkSymmetry()

print "\nIterative Solver (1)=>"
u_i=mypde.getSolution()

print "\nDirect Solver (1)=>"
mypde.setSolverMethod(DIRECT)
u_d=mypde.getSolution()


print "\n***************************************************************"
error=u_ex-u_d
print "norm of the error for direct solver is   : ",error.Lsup()/norm_u_ex
error=u_ex-u_i
print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex
print "***************************************************************"

# get handles to nodes and elements 2

print "\nGet handles to nodes and elements(2)=>"
print "--------------------------------------"

print "x=n.getX():"
x=n.getX()

print "msk=x[0].whereZero()+(x[0]-1.).whereZero()"
msk=x[0].whereZero()+(x[0]-1.).whereZero()

print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
mypde=LinearPDE(mydomain)
mypde.setDebugOn()
mypde.setValue(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)

print "mypde.checkSymmetry()"
print mypde.checkSymmetry()

# generate a test solution 2

print "\nGenerate a test solution (2)"
print "----------------------------"

print "\nDirect Solver (2)=>"

# mypde.setSymmetryOn() : is not woking yet!
mypde.setSolverMethod(mypde.DIRECT)
u_d=mypde.getSolution()

print "\nIterative Solver (2)=>"

mypde.setSymmetryOn() 
mypde.setSolverMethod(DEFAULT_METHOD)
u_i=mypde.getSolution()

print "\n******************************************************************"
error=u_ex-u_d
print "norm of the error for direct solver is   : ",error.Lsup()/norm_u_ex
error=u_ex-u_i
print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex
print "******************************************************************"

print "\n-----"
print "Done."
print "-----"

grad(u_d).saveDX("t.dx")
1	jgs	82	# $Id$
2
3			import sys
4			import os
5			import unittest
6
7	jgs	102	from esys.escript import *
8			from esys.linearPDEs import *
9			from esys import finley
10	jgs	82
11			print "\nSimpleSolve.py"
12			print "--------------"
13
14			alpha=0.01
15
16			# generate mesh
17
18	jgs	102	# print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
19			# mydomain=finley.Rectangle(1,1)
20	jgs	82
21	jgs	102	# print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
22	jgs	113	mydomain=finley.Rectangle(4,4,1)
23	jgs	82
24			print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
25	jgs	113	# mydomain=finley.Rectangle(151,151,1)
26	jgs	82
27			print "\nSetup domain and functions"
28			print "--------------------------"
29
30			print "e=Function(mydomain):"
31			e=Function(mydomain)
32
33			print "n=ContinuousFunction(mydomain):"
34			n=ContinuousFunction(mydomain)
35
36			# get handles to nodes and elements 1
37
38			print "\nGet handles to nodes and elements(1)=>"
39			print "--------------------------------------"
40
41			print "u_ex=Scalar(1,n,True):"
42			u_ex=Scalar(1,n,True)
43
44			print "x=e.getX():"
45			x=e.getX()
46
47			print "norm_u_ex=u_ex.Lsup():"
48			norm_u_ex=u_ex.Lsup()
49
50	jgs	108	print "\nGenerate a test solution (1)"
51			print "----------------------------"
52
53	jgs	102	print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
54	jgs	104	mypde=LinearPDE(mydomain)
55	jgs	102	mypde.setDebugOn()
56	jgs	104	mypde.setValue(A=[[1.,0.8],[0.4,1.]],D=alpha,Y=alpha)
57	jgs	82
58	jgs	108	print "mypde.checkSymmetry()"
59			print mypde.checkSymmetry()
60	jgs	82
61	jgs	102	print "\nIterative Solver (1)=>"
62			u_i=mypde.getSolution()
63
64	jgs	100	print "\nDirect Solver (1)=>"
65	jgs	102	mypde.setSolverMethod(DIRECT)
66			u_d=mypde.getSolution()
67	jgs	97
68	jgs	82
69			print "\n***************************************************************"
70			error=u_ex-u_d
71			print "norm of the error for direct solver is : ",error.Lsup()/norm_u_ex
72			error=u_ex-u_i
73			print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex
74			print "***************************************************************"
75
76			# get handles to nodes and elements 2
77
78			print "\nGet handles to nodes and elements(2)=>"
79			print "--------------------------------------"
80
81			print "x=n.getX():"
82			x=n.getX()
83
84			print "msk=x[0].whereZero()+(x[0]-1.).whereZero()"
85			msk=x[0].whereZero()+(x[0]-1.).whereZero()
86
87	jgs	102	print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
88	jgs	104	mypde=LinearPDE(mydomain)
89	jgs	108	mypde.setDebugOn()
90	jgs	104	mypde.setValue(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)
91	jgs	82
92	jgs	108	print "mypde.checkSymmetry()"
93			print mypde.checkSymmetry()
94
95	jgs	82	# generate a test solution 2
96
97			print "\nGenerate a test solution (2)"
98			print "----------------------------"
99
100			print "\nDirect Solver (2)=>"
101
102	jgs	102	# mypde.setSymmetryOn() : is not woking yet!
103	jgs	104	mypde.setSolverMethod(mypde.DIRECT)
104	jgs	102	u_d=mypde.getSolution()
105	jgs	82
106			print "\nIterative Solver (2)=>"
107
108	jgs	108	mypde.setSymmetryOn()
109	jgs	102	mypde.setSolverMethod(DEFAULT_METHOD)
110			u_i=mypde.getSolution()
111	jgs	82
112			print "\n******************************************************************"
113			error=u_ex-u_d
114			print "norm of the error for direct solver is : ",error.Lsup()/norm_u_ex
115			error=u_ex-u_i
116			print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex
117			print "******************************************************************"
118
119			print "\n-----"
120			print "Done."
121			print "-----"
122	jgs	113
123			grad(u_d).saveDX("t.dx")
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision