| 1 |
# $Id$ |
| 2 |
|
| 3 |
__copyright__=""" Copyright (c) 2006 by ACcESS MNRF |
| 4 |
http://www.access.edu.au |
| 5 |
Primary Business: Queensland, Australia""" |
| 6 |
__license__="""Licensed under the Open Software License version 3.0 |
| 7 |
http://www.opensource.org/licenses/osl-3.0.php""" |
| 8 |
import sys |
| 9 |
import os |
| 10 |
import unittest |
| 11 |
import time |
| 12 |
|
| 13 |
from esys.escript import * |
| 14 |
from esys.escript.linearPDEs import * |
| 15 |
from esys import finley |
| 16 |
|
| 17 |
starttime = time.clock() |
| 18 |
|
| 19 |
print "\nSimpleSolve.py" |
| 20 |
print "--------------" |
| 21 |
|
| 22 |
alpha=0.7 |
| 23 |
error_tol=1.e-5 |
| 24 |
|
| 25 |
# generate mesh |
| 26 |
|
| 27 |
# print "\nGenerate mesh: finley.Rectangle(9,12,1)=>" |
| 28 |
# mydomain=finley.Rectangle(140,140) |
| 29 |
|
| 30 |
# print "\nGenerate mesh: finley.Rectangle(4,4,1)=>" |
| 31 |
# mydomain=finley.Rectangle(50,50,1) |
| 32 |
mydomain=finley.Rectangle(500,500,1) |
| 33 |
# mydomain=finley.Rectangle(150,150,1) |
| 34 |
|
| 35 |
print "\nGenerate mesh: finley.Rectangle(151,151,1)=>" |
| 36 |
# mydomain=finley.Rectangle(151,151,1) |
| 37 |
# mydomain=finley.Rectangle(128,128,1) |
| 38 |
|
| 39 |
# set the direct solver switch |
| 40 |
DIRECT=LinearPDE.DIRECT |
| 41 |
# DIRECT=LinearPDE.ITERATIVE # this will switch of the DIRECT solver call to avoid external library calls which may not be available everywhere |
| 42 |
|
| 43 |
print "\nSetup domain and functions" |
| 44 |
print "--------------------------" |
| 45 |
|
| 46 |
print "e=Function(mydomain):" |
| 47 |
e=Function(mydomain) |
| 48 |
|
| 49 |
print "n=ContinuousFunction(mydomain):" |
| 50 |
n=ContinuousFunction(mydomain) |
| 51 |
|
| 52 |
# get handles to nodes and elements 1 |
| 53 |
|
| 54 |
print "\nGet handles to nodes and elements(1)=>" |
| 55 |
print "--------------------------------------" |
| 56 |
|
| 57 |
print "u_ex=Scalar(1,n,True):" |
| 58 |
u_ex=Scalar(1,n,True) |
| 59 |
|
| 60 |
print "x=e.getX():" |
| 61 |
x=e.getX() |
| 62 |
|
| 63 |
print "norm_u_ex=Lsup(u_ex):" |
| 64 |
norm_u_ex=Lsup(u_ex) |
| 65 |
|
| 66 |
print "\nGenerate a test solution (1)" |
| 67 |
print "----------------------------" |
| 68 |
|
| 69 |
print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)" |
| 70 |
mypde=LinearPDE(mydomain) |
| 71 |
mypde.setDebugOn() |
| 72 |
mypde.setValue(A=[[1.,-0.001],[-0.001,1.]],D=alpha,Y=alpha) |
| 73 |
|
| 74 |
print "mypde.checkSymmetry()" |
| 75 |
print mypde.checkSymmetry() |
| 76 |
|
| 77 |
print "\nIterative Solver (1)=>" |
| 78 |
# mypde.setSolverMethod(mypde.PRES20,preconditioner=mypde.ILU0) |
| 79 |
mypde.setSolverMethod(mypde.BICGSTAB,preconditioner=mypde.JACOBI) |
| 80 |
u_i=mypde.getSolution(verbose=True,iter_max=3000) |
| 81 |
|
| 82 |
print "\nDirect Solver (1)=>" |
| 83 |
mypde.setSolverMethod(DIRECT) |
| 84 |
u_d=mypde.getSolution(verbose=True) |
| 85 |
|
| 86 |
print "\n***************************************************************" |
| 87 |
error=u_ex-u_d |
| 88 |
error_norm=Lsup(error)/norm_u_ex |
| 89 |
print "norm of the error for direct solver is : ",error_norm |
| 90 |
if error_norm > error_tol: |
| 91 |
print "### error norm exceeded maximum tolerance ###" |
| 92 |
sys.exit(1) |
| 93 |
error=u_ex-u_i |
| 94 |
error_norm=Lsup(error)/norm_u_ex |
| 95 |
print "norm of the error for iterative solver is: ",error_norm |
| 96 |
if error_norm > error_tol: |
| 97 |
print "### error norm exceeded maximum tolerance ###" |
| 98 |
sys.exit(1) |
| 99 |
print "***************************************************************" |
| 100 |
del mypde |
| 101 |
print "***************************************************************" |
| 102 |
|
| 103 |
|
| 104 |
# get handles to nodes and elements 2 |
| 105 |
|
| 106 |
print "\nGet handles to nodes and elements(2)=>" |
| 107 |
print "--------------------------------------" |
| 108 |
|
| 109 |
print "x=n.getX():" |
| 110 |
x=n.getX() |
| 111 |
|
| 112 |
print " msk=whereZero(x[0])+whereZero(x[0]-1.)" |
| 113 |
msk=whereZero(x[0])+whereZero(x[0]-1.) |
| 114 |
|
| 115 |
print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)" |
| 116 |
mypde=LinearPDE(mydomain) |
| 117 |
mypde.setDebugOn() |
| 118 |
mypde.setValue(A=[[1.,0.0],[0.0,1.]],q=msk,r=u_ex) |
| 119 |
|
| 120 |
print "mypde.checkSymmetry()" |
| 121 |
print mypde.checkSymmetry() |
| 122 |
|
| 123 |
# generate a test solution 2 |
| 124 |
|
| 125 |
print "\nGenerate a test solution (2)" |
| 126 |
print "----------------------------" |
| 127 |
|
| 128 |
print "\nDirect Solver (2)=>" |
| 129 |
|
| 130 |
mypde.setSymmetryOn() |
| 131 |
mypde.setTolerance(1.e-13) |
| 132 |
|
| 133 |
# mypde.setSymmetryOn() : is not woking yet! |
| 134 |
mypde.setSolverMethod(DIRECT) |
| 135 |
u_d=mypde.getSolution(verbose=True) |
| 136 |
|
| 137 |
print "\nIterative Solver (2)=>" |
| 138 |
|
| 139 |
mypde.setSolverMethod(mypde.PCG) |
| 140 |
u_i=mypde.getSolution(verbose=True,iter_max=3000) |
| 141 |
|
| 142 |
print "\n******************************************************************" |
| 143 |
error=u_ex-u_d |
| 144 |
error_norm=Lsup(error)/norm_u_ex |
| 145 |
print "norm of the error for direct solver is : ",error_norm |
| 146 |
if error_norm > error_tol: |
| 147 |
print "### error norm exceeded maximum tolerance ###" |
| 148 |
sys.exit(1) |
| 149 |
error=u_ex-u_i |
| 150 |
error_norm=Lsup(error)/norm_u_ex |
| 151 |
print "norm of the error for iterative solver is: ",error_norm |
| 152 |
if error_norm > error_tol: |
| 153 |
print "### error norm exceeded maximum tolerance ###" |
| 154 |
sys.exit(1) |
| 155 |
print "******************************************************************" |
| 156 |
|
| 157 |
print "\n-----" |
| 158 |
print "Done." |
| 159 |
print "-----" |
| 160 |
|
| 161 |
stoptime = time.clock() |
| 162 |
elapsed = stoptime - starttime |
| 163 |
print "\nElapsed time: ", elapsed, "\n" |
| 164 |
|
| 165 |
sys.exit(0) |
Parent Directory
| 
Revision Log