| 1 |
# $Id$ |
| 2 |
|
| 3 |
import sys |
| 4 |
import os |
| 5 |
import unittest |
| 6 |
|
| 7 |
from esys.escript import * |
| 8 |
from esys.linearPDEs import * |
| 9 |
from esys import finley |
| 10 |
|
| 11 |
print "\nSimpleSolve.py" |
| 12 |
print "--------------" |
| 13 |
|
| 14 |
alpha=0.01 |
| 15 |
|
| 16 |
# generate mesh |
| 17 |
|
| 18 |
# print "\nGenerate mesh: finley.Rectangle(9,12,1)=>" |
| 19 |
# mydomain=finley.Rectangle(1,1) |
| 20 |
|
| 21 |
# print "\nGenerate mesh: finley.Rectangle(4,4,1)=>" |
| 22 |
# mydomain=finley.Rectangle(4,4,1) |
| 23 |
|
| 24 |
print "\nGenerate mesh: finley.Rectangle(151,151,1)=>" |
| 25 |
mydomain=finley.Rectangle(151,151,1) |
| 26 |
|
| 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 |
print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)" |
| 51 |
mypde=LinearPDE(A=[[1.,0.8],[0.4,1.]],D=alpha,Y=alpha,domain=mydomain) |
| 52 |
mypde.setDebugOn() |
| 53 |
#mypde=LinearPDE(A=[[1.,0.],[0.,1.]],D=alpha,Y=alpha,domain=mydomain) |
| 54 |
mypde.getOperator().saveMM("t.msh") |
| 55 |
|
| 56 |
# generate a test solution 1 |
| 57 |
|
| 58 |
print "\nGenerate a test solution (1)" |
| 59 |
print "----------------------------" |
| 60 |
|
| 61 |
print "\nIterative Solver (1)=>" |
| 62 |
|
| 63 |
u_i=mypde.getSolution() |
| 64 |
|
| 65 |
print "\nDirect Solver (1)=>" |
| 66 |
|
| 67 |
mypde.setSolverMethod(DIRECT) |
| 68 |
u_d=mypde.getSolution() |
| 69 |
|
| 70 |
|
| 71 |
print "\n***************************************************************" |
| 72 |
error=u_ex-u_d |
| 73 |
print "norm of the error for direct solver is : ",error.Lsup()/norm_u_ex |
| 74 |
error=u_ex-u_i |
| 75 |
print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex |
| 76 |
print "***************************************************************" |
| 77 |
|
| 78 |
# get handles to nodes and elements 2 |
| 79 |
|
| 80 |
print "\nGet handles to nodes and elements(2)=>" |
| 81 |
print "--------------------------------------" |
| 82 |
|
| 83 |
print "x=n.getX():" |
| 84 |
x=n.getX() |
| 85 |
|
| 86 |
print "msk=x[0].whereZero()+(x[0]-1.).whereZero()" |
| 87 |
msk=x[0].whereZero()+(x[0]-1.).whereZero() |
| 88 |
|
| 89 |
print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)" |
| 90 |
mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex) |
| 91 |
mypde.setDebugOn() |
| 92 |
|
| 93 |
# generate a test solution 2 |
| 94 |
|
| 95 |
print "\nGenerate a test solution (2)" |
| 96 |
print "----------------------------" |
| 97 |
|
| 98 |
print "\nDirect Solver (2)=>" |
| 99 |
|
| 100 |
# mypde.setSymmetryOn() : is not woking yet! |
| 101 |
mypde.setSolverMethod(DIRECT) |
| 102 |
u_d=mypde.getSolution() |
| 103 |
|
| 104 |
print "\nIterative Solver (2)=>" |
| 105 |
|
| 106 |
#mypde.setSymmetryOn() |
| 107 |
mypde.setSolverMethod(DEFAULT_METHOD) |
| 108 |
u_i=mypde.getSolution() |
| 109 |
|
| 110 |
print "\n******************************************************************" |
| 111 |
error=u_ex-u_d |
| 112 |
print "norm of the error for direct solver is : ",error.Lsup()/norm_u_ex |
| 113 |
error=u_ex-u_i |
| 114 |
print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex |
| 115 |
print "******************************************************************" |
| 116 |
|
| 117 |
print "\n-----" |
| 118 |
print "Done." |
| 119 |
print "-----" |
Parent Directory
| 
Revision Log