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

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

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

revision 115 by jgs, Fri Mar 4 07:12:47 2005 UTC revision 117 by jgs, Fri Apr 1 05:48:57 2005 UTC
# Line 17  Line 17 
17  from esys.escript import *  from esys.escript import *
18  from esys.linearPDEs import *  from esys.linearPDEs import *
19  import esys.finley as pdelib  import esys.finley as pdelib
20    from time import time
21    
22  from numarray import *  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,prec):  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 74  def solveVector(numDim, totalNumElem, le Line 81  def solveVector(numDim, totalNumElem, le
81      # Build the pdelib System Matrix and RHS      # Build the pdelib System Matrix and RHS
82      mypde=LinearPDE(mesh)      mypde=LinearPDE(mesh)
83      mypde.setValue(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    
86      # Solve for Approximate Solution      # Solve for Approximate Solution
87      u_approx = mypde.getSolution(preconditioner=prec,iter_max=10000)      tm=time()
88        u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000)
89        tm=time()-tm
90    
91      # Report Results      # Report Results
92      error=Lsup(u - u_approx)/Lsup(u)      error=Lsup(u - u_approx)/Lsup(u)
93      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))  
94        
95      return error      return error
96    
# Line 90  def solveVector(numDim, totalNumElem, le Line 98  def solveVector(numDim, totalNumElem, le
98    
99  def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):  def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
100    
101        if prec=="":
102            prec_id=0
103        else:
104           prec_id=eval("LinearPDE.%s"%prec)
105        solver_method_id=eval("LinearPDE.%s"%solver_method)
106      print "Scalar solver:"      print "Scalar solver:"
107      recDim=array([len_x0,1.,1.])      recDim=array([len_x0,1.,1.])
108      # Define Computational Domain      # Define Computational Domain
# Line 128  def solveScalar(numDim, totalNumElem, le Line 141  def solveScalar(numDim, totalNumElem, le
141      # Build the pdelib System Matrix and RHS      # Build the pdelib System Matrix and RHS
142      mypde=LinearPDE(mesh)      mypde=LinearPDE(mesh)
143      mypde.setValue(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)
144      mypde.setSolverMethod(solver_method)      mypde.setSolverMethod(solver_method_id)
145    
146      # Solve for Approximate Solution      # Solve for Approximate Solution
147      u_approx = mypde.getSolution(preconditioner=prec,iter_max=10000)      tm=time()
148        u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000)
149        tm=time()-tm
150    
151      # Report Results      # Report Results
152      error=Lsup(u - u_approx)/Lsup(u)      error=Lsup(u - u_approx)/Lsup(u)
153      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))  
154        
155      return error      return error
156    
# Line 146  def solveScalar(numDim, totalNumElem, le Line 160  def solveScalar(numDim, totalNumElem, le
160  print "Test is started:"  print "Test is started:"
161  print "----------------"  print "----------------"
162  error=0.  error=0.
163  # for numDim in [2,3]:  for numDim in [2,3]:
164  for numDim in [3]:     for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
165     for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400]:        for problem in [solveScalar,solveVector]:
166        # for problem in [solveScalar,solveVector]:           if totalNumElem*2**numDim*numDim< 200000: error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
167        for problem in [solveScalar]:           for solver_method in [ "PCG" ]:
168           # for solver_method in [ LinearPDE.PRES20, LinearPDE.PCG, LinearPDE.DIRECT, LinearPDE.BICGSTAB]:              for prec in [ "JACOBI", "ILU0" ]:
          for solver_method in [ LinearPDE.PCG ]:  
             # for prec in [ LinearPDE.JACOBI, LinearPDE.ILU0 ]:  
             for prec in [ LinearPDE.ILU0 ]:  
169                 error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error])                 error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error])
170  print "----------------"  print "----------------"
171  print "maximum error over all tests is ",error  print "maximum error over all tests is ",error
Legend:
Removed from v.115  
changed lines
  Added in v.117
  ViewVC Help
Powered by ViewVC 1.1.26