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

Revision 122 - (show annotations)
Thu Jun 9 05:38:05 2005 UTC (17 years, 9 months ago) by jgs
File MIME type: text/x-python
File size: 6078 byte(s)
```Merge of development branch back to main trunk on 2005-06-09

```
 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 from time import time 21 22 from numarray import * 23 24 # these values are currently fixed: 25 len_x0=1. 26 alpha=10. 27 tm=0 28 29 #############################################################################################################3 30 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:" 38 recDim=array([len_x0,1.,1.]) 39 # Define Computational Domain 40 numElem=int((totalNumElem/(len_x0*1.))**(1./numDim)) 41 elemDim = array([int(len_x0*numElem), numElem, numElem],Int) 42 43 # Set Mesh 44 if (numDim == 2): 45 mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \ 46 l0 = len_x0, l1 = 1.) 47 totElem=elemDim[0]*elemDim[1] 48 elif (numDim == 3): 49 mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \ 50 l0 = len_x0, l1 = 1., l2 = 1.) 51 totElem=elemDim[0]*elemDim[1]*elemDim[2] 52 53 print " length of domain: ",recDim[:numDim] 54 print " requested elements: ",totalNumElem 55 print " num elements: ",totElem 56 # Set Mesh Descriptors 57 meshDim = mesh.getDim() 58 contfunc = ContinuousFunction(mesh) 59 func = Function(mesh) 60 x = contfunc.getX() 61 62 # Set Boundary Mask / pdelib Template "q" Parameter Vector 63 bndryMask = Vector(value = 0, what = contfunc) 64 for i in range(meshDim): 65 bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) \ 66 * ones((numDim,)) 67 68 # Set True Solution / pdelib Template "r" Parameter Vector 69 u = Vector(value = 0, what = contfunc) 70 for i in range(meshDim): 71 for j in range(meshDim - 1): 72 u[i] += x[(i + j + 1) % meshDim]**2 73 # Set pdelib Template "A" Parameter Tensor 74 A = Tensor4(value = 0, what = func) 75 for i in range(meshDim): 76 for j in range(meshDim): 77 A[i,i,j,j] += 1. 78 A[i,j,j,i] += alpha 79 A[i,j,i,j] += alpha 80 81 # Build the pdelib System Matrix and RHS 82 mypde=LinearPDE(mesh) 83 mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u) 84 mypde.setSolverMethod(solver_method_id) 85 86 # Solve for Approximate Solution 87 tm=time() 88 u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000) 89 tm=time()-tm 90 91 # Report Results 92 error=Lsup(u - u_approx)/Lsup(u) 93 print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm) 94 95 return error 96 97 ################################################################################################################# 98 99 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:" 107 recDim=array([len_x0,1.,1.]) 108 # Define Computational Domain 109 numElem=int((totalNumElem/(len_x0*1.))**(1./numDim)) 110 elemDim = array([int(len_x0*numElem), numElem, numElem],Int) 111 # Set Mesh 112 if (numDim == 2): 113 mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \ 114 l0 = len_x0, l1 = 1.) 115 totElem=elemDim[0]*elemDim[1] 116 elif (numDim == 3): 117 mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \ 118 l0 = len_x0, l1 = 1., l2 = 1.) 119 totElem=elemDim[0]*elemDim[1]*elemDim[2] 120 121 print " length of domain: ",recDim[:numDim] 122 print " requested elements: ",totalNumElem 123 print " num elements: ",totElem 124 125 # Set Mesh Descriptors 126 meshDim = mesh.getDim() 127 contfunc = ContinuousFunction(mesh) 128 func = Function(mesh) 129 x = contfunc.getX() 130 131 # Set Boundary Mask / pdelib Template "q" Parameter Vector 132 bndryMask = Scalar(value = 0, what = contfunc) 133 for i in range(meshDim): 134 bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) * 1.0 135 136 # Set True Solution / pdelib Template "r" Parameter Vector 137 u = Scalar(value = 0, what = contfunc) 138 for j in range(meshDim): 139 u += x[j] * x[j] 140 141 # Build the pdelib System Matrix and RHS 142 mypde=LinearPDE(mesh) 143 mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u) 144 mypde.setSolverMethod(solver_method_id) 145 146 # Solve for Approximate Solution 147 tm=time() 148 u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000) 149 tm=time()-tm 150 151 # Report Results 152 error=Lsup(u - u_approx)/Lsup(u) 153 print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm) 154 155 return error 156 157 ####################################################################################### 158 159 160 print "Test is started:" 161 print "----------------" 162 error=0. 163 for numDim in [2, 3]: 164 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]: 165 for problem in [solveScalar,solveVector]: 166 error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error]) 167 #if totalNumElem*2**numDim*numDim< 200000: error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error]) 168 # for solver_method in [ "PCG" ]: 169 # for prec in [ "JACOBI", "ILU0" ]: 170 # error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error]) 171 print "----------------" 172 print "maximum error over all tests is ",error 173 print "----------------"

## Properties

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