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

Revision 115 - (show annotations)
Fri Mar 4 07:12:47 2005 UTC (15 years, 5 months ago) by jgs
File MIME type: text/x-python
File size: 5683 byte(s)
```*** empty log message ***

```
 1 # \$Id\$ 2 3 """General test environment to test the solvers for scalar and vector equations 4 5 test parameters are 6 7 numDim = spatial dimension 8 totalNumElem = number of func in each direction 9 problem = solveScalar,solveVector 10 11 solver_method = true/false 12 len_x0 = length of the domain in x0 direction (number of func in x0 is round(totalNumElem*len_x0) ) 13 alpha = a parameter of the PDE (not well defined yet) 14 15 """ 16 17 from esys.escript import * 18 from esys.linearPDEs import * 19 import esys.finley as pdelib 20 21 from numarray import * 22 23 # these values are currently fixed: 24 len_x0=1. 25 alpha=0.1 26 27 #############################################################################################################3 28 def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec): 29 30 print "Vector solver:" 31 recDim=array([len_x0,1.,1.]) 32 # Define Computational Domain 33 numElem=int((totalNumElem/(len_x0*1.))**(1./numDim)) 34 elemDim = array([int(len_x0*numElem), numElem, numElem],Int) 35 36 # Set Mesh 37 if (numDim == 2): 38 mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \ 39 l0 = len_x0, l1 = 1.) 40 totElem=elemDim[0]*elemDim[1] 41 elif (numDim == 3): 42 mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \ 43 l0 = len_x0, l1 = 1., l2 = 1.) 44 totElem=elemDim[0]*elemDim[1]*elemDim[2] 45 46 print " length of domain: ",recDim[:numDim] 47 print " requested elements: ",totalNumElem 48 print " num elements: ",totElem 49 # Set Mesh Descriptors 50 meshDim = mesh.getDim() 51 contfunc = ContinuousFunction(mesh) 52 func = Function(mesh) 53 x = contfunc.getX() 54 55 # Set Boundary Mask / pdelib Template "q" Parameter Vector 56 bndryMask = Vector(value = 0, what = contfunc) 57 for i in range(meshDim): 58 bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) \ 59 * ones((numDim,)) 60 61 # Set True Solution / pdelib Template "r" Parameter Vector 62 u = Vector(value = 0, what = contfunc) 63 for i in range(meshDim): 64 for j in range(meshDim - 1): 65 u[i] += x[(i + j + 1) % meshDim]**2 66 # Set pdelib Template "A" Parameter Tensor 67 A = Tensor4(value = 0, what = func) 68 for i in range(meshDim): 69 for j in range(meshDim): 70 A[i,i,j,j] += 1. 71 A[i,j,j,i] += alpha 72 A[i,j,i,j] += alpha 73 74 # Build the pdelib System Matrix and RHS 75 mypde=LinearPDE(mesh) 76 mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u) 77 mypde.setSolverMethod(solver_method) 78 79 # Solve for Approximate Solution 80 u_approx = mypde.getSolution(preconditioner=prec,iter_max=10000) 81 82 # Report Results 83 error=Lsup(u - u_approx)/Lsup(u) 84 print " error L^sup Norm : ", error 85 print " residual L^sup Norm : ", Lsup(mypde.getResidual(u_approx)) 86 87 return error 88 89 ################################################################################################################# 90 91 def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec): 92 93 print "Scalar solver:" 94 recDim=array([len_x0,1.,1.]) 95 # Define Computational Domain 96 numElem=int((totalNumElem/(len_x0*1.))**(1./numDim)) 97 elemDim = array([int(len_x0*numElem), numElem, numElem],Int) 98 # Set Mesh 99 if (numDim == 2): 100 mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \ 101 l0 = len_x0, l1 = 1.) 102 totElem=elemDim[0]*elemDim[1] 103 elif (numDim == 3): 104 mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \ 105 l0 = len_x0, l1 = 1., l2 = 1.) 106 totElem=elemDim[0]*elemDim[1]*elemDim[2] 107 108 print " length of domain: ",recDim[:numDim] 109 print " requested elements: ",totalNumElem 110 print " num elements: ",totElem 111 112 # Set Mesh Descriptors 113 meshDim = mesh.getDim() 114 contfunc = ContinuousFunction(mesh) 115 func = Function(mesh) 116 x = contfunc.getX() 117 118 # Set Boundary Mask / pdelib Template "q" Parameter Vector 119 bndryMask = Scalar(value = 0, what = contfunc) 120 for i in range(meshDim): 121 bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) * 1.0 122 123 # Set True Solution / pdelib Template "r" Parameter Vector 124 u = Scalar(value = 0, what = contfunc) 125 for j in range(meshDim): 126 u += x[j] * x[j] 127 128 # Build the pdelib System Matrix and RHS 129 mypde=LinearPDE(mesh) 130 mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u) 131 mypde.setSolverMethod(solver_method) 132 133 # Solve for Approximate Solution 134 u_approx = mypde.getSolution(preconditioner=prec,iter_max=10000) 135 136 # Report Results 137 error=Lsup(u - u_approx)/Lsup(u) 138 print " error L^sup Norm : ", error 139 print " residual L^sup Norm : ", Lsup(mypde.getResidual(u_approx)) 140 141 return error 142 143 ####################################################################################### 144 145 146 print "Test is started:" 147 print "----------------" 148 error=0. 149 # for numDim in [2,3]: 150 for numDim in [3]: 151 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400]: 152 # for problem in [solveScalar,solveVector]: 153 for problem in [solveScalar]: 154 # for solver_method in [ LinearPDE.PRES20, LinearPDE.PCG, LinearPDE.DIRECT, LinearPDE.BICGSTAB]: 155 for solver_method in [ LinearPDE.PCG ]: 156 # for prec in [ LinearPDE.JACOBI, LinearPDE.ILU0 ]: 157 for prec in [ LinearPDE.ILU0 ]: 158 error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error]) 159 print "----------------" 160 print "maximum error over all tests is ",error 161 print "----------------"

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision