/[escript]/trunk/esys2/finley/test/python/SolveTest.py
ViewVC logotype

Diff of /trunk/esys2/finley/test/python/SolveTest.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 102 by jgs, Wed Dec 15 07:08:39 2004 UTC revision 153 by jgs, Tue Oct 25 01:51:20 2005 UTC
# Line 14  Line 14 
14    
15  """  """
16    
17  import sys  from esys.escript import *
18  import os  from esys.escript.linearPDEs import *
19                                                                                                                                                        import esys.finley as pdelib
20  esys_root=os.getenv('ESYS_ROOT')  from time import time
 sys.path.append(esys_root+'/finley/lib')  
 sys.path.append(esys_root+'/escript/lib')  
 sys.path.append(esys_root+'/escript/py_src')  
                                                                                                                                                       
 from escript import *  
 from util import *  
 from linearPDEs import *  
 from numarray import *  
 import finley as pdelib  
21    
22    from numarray import *
23    
24  # these values are currently fixed:  # these values are currently fixed:
25  len_x0=1.  len_x0=1.
26  alpha=0.1  alpha=10.
27    tm=0
28    
29  #############################################################################################################3  #############################################################################################################3
30  def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method):  def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
31        
32        if prec=="":
33            prec_id=0
34        else:
35           prec_id=eval("LinearPDE.%s"%prec)
36        solver_method_id=eval("LinearPDE.%s"%solver_method)
37      print "Vector solver:"      print "Vector solver:"
38      recDim=array([len_x0,1.,1.])      recDim=array([len_x0,1.,1.])
39      # Define Computational Domain      # Define Computational Domain
# Line 81  def solveVector(numDim, totalNumElem, le Line 79  def solveVector(numDim, totalNumElem, le
79          A[i,j,i,j] += alpha          A[i,j,i,j] += alpha
80    
81      # Build the pdelib System Matrix and RHS      # Build the pdelib System Matrix and RHS
82      mypde=LinearPDE(domain=mesh, \      mypde=LinearPDE(mesh)
83                       A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)      mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)
84      mypde.setSolverMethod(solver_method)      mypde.setSolverMethod(solver_method_id)
85        # mypde.getOperator().saveMM("g.mm")
86    
87      # Solve for Approximate Solution      # Solve for Approximate Solution
88      u_approx = mypde.getSolution(iter_max=10000)      tm=time()
89        u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
90        tm=time()-tm
91    
92      # Report Results      # Report Results
93      error=Lsup(u - u_approx)/Lsup(u)      error=Lsup(u - u_approx)/Lsup(u)
94      print "    error L^sup Norm : ", error      print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
     print "    residual L^sup Norm : ", Lsup(mypde.getResidual(u_approx))  
95        
96      return error      return error
97    
98  #################################################################################################################  #################################################################################################################
99    
100  def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method):  def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
101    
102        if prec=="":
103            prec_id=0
104        else:
105           prec_id=eval("LinearPDE.%s"%prec)
106        solver_method_id=eval("LinearPDE.%s"%solver_method)
107      print "Scalar solver:"      print "Scalar solver:"
108      recDim=array([len_x0,1.,1.])      recDim=array([len_x0,1.,1.])
109      # Define Computational Domain      # Define Computational Domain
# Line 135  def solveScalar(numDim, totalNumElem, le Line 140  def solveScalar(numDim, totalNumElem, le
140          u += x[j] * x[j]          u += x[j] * x[j]
141    
142      # Build the pdelib System Matrix and RHS      # Build the pdelib System Matrix and RHS
143      mypde=LinearPDE(domain=mesh, \      mypde=LinearPDE(mesh)
144                      A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)      mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
145      mypde.setSolverMethod(solver_method)      mypde.setSolverMethod(solver_method_id)
146    
147      # Solve for Approximate Solution      # Solve for Approximate Solution
148      u_approx = mypde.getSolution(iter_max=10000)      tm=time()
149        u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
150        tm=time()-tm
151    
152      # Report Results      # Report Results
153      error=Lsup(u - u_approx)/Lsup(u)      error=Lsup(u - u_approx)/Lsup(u)
154      print "    error L^sup Norm : ", error      print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
     print "    residual L^sup Norm : ", Lsup(mypde.getResidual(u_approx))  
155        
156      return error      return error
157    
# Line 155  def solveScalar(numDim, totalNumElem, le Line 161  def solveScalar(numDim, totalNumElem, le
161  print "Test is started:"  print "Test is started:"
162  print "----------------"  print "----------------"
163  error=0.  error=0.
164  for numDim in [2,3]:  for numDim in [2, 3]:
165     for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400]:     # for totalNumElem in [51200]:
166       for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
167        for problem in [solveScalar,solveVector]:        for problem in [solveScalar,solveVector]:
168           # for solver_method in [ LinearPDE.PRES20, LinearPDE.PCG, LinearPDE.DIRECT, LinearPDE.BICGSTAB]:        #for problem in [solveVector]:
169           for solver_method in [ LinearPDE.PRES20, LinearPDE.PCG, LinearPDE.BICGSTAB]:           error=max([problem(numDim, totalNumElem, len_x0, alpha,"PCG",""),error])
170              error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method),error])           error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
171             #if totalNumElem*2**numDim*numDim< 200000: error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
172             # for solver_method in [ "PCG" ]:
173             #    for prec in [ "JACOBI", "ILU0" ]:
174             #       error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error])
175  print "----------------"  print "----------------"
176  print "maximum error over all tests is ",error  print "maximum error over all tests is ",error
177  print "----------------"  print "----------------"

Legend:
Removed from v.102  
changed lines
  Added in v.153

  ViewVC Help
Powered by ViewVC 1.1.26