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

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

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

revision 3096 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 3097 by gross, Fri Aug 20 04:59:12 2010 UTC
# Line 976  class SimpleSolve_Brick_Order2_SystemPDE Line 976  class SimpleSolve_Brick_Order2_SystemPDE
976          error=Lsup(u-u_ex)          error=Lsup(u-u_ex)
977          self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)          self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
978    
979    class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
980         def test_solve(self):
981            domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
982            x=Solution(domain).getX()
983            # --- set exact solution ----
984            u_ex=Scalar(0,Solution(domain))
985            u_ex=1.+2.*x[0]+3.*x[1]
986            # --- set exact gradient -----------
987            g_ex=Data(0.,(2,),Solution(domain))
988            g_ex[0]=2.
989            g_ex[1]=3.
990            # -------- test gradient --------------------------------
991            g=grad(u_ex)
992            self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
993            # -------- set-up PDE -----------------------------------
994            pde=LinearPDE(domain,numEquations=1)
995            mask=whereZero(x[0])
996            pde.setValue(r=u_ex,q=mask)
997            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
998            # -------- get the solution ---------------------------
999            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1000            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1001            pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1002            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1003            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1004            u=pde.getSolution()
1005            # -------- test the solution ---------------------------
1006            error=Lsup(u-u_ex)
1007            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1008    
1009    class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1010         def test_solve(self):
1011            domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
1012            x=Solution(domain).getX()
1013            # --- set exact solution ----
1014            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1015            # --- set exact gradient -----------
1016            g_ex=Data(0.,(2,),Solution(domain))
1017            g_ex[0]=2.+8.*x[0]+5.*x[1]
1018            g_ex[1]=3.+5.*x[0]+12.*x[1]
1019            # -------- test gradient --------------------------------
1020            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1021            # -------- set-up PDE -----------------------------------
1022            pde=LinearPDE(domain,numEquations=1)
1023            mask=whereZero(x[0])
1024            pde.setValue(r=u_ex,q=mask)
1025            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
1026            # -------- get the solution ---------------------------
1027            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1028            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1029            pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1030            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1031            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1032            u=pde.getSolution()
1033            # -------- test the solution ---------------------------
1034            error=Lsup(u-u_ex)
1035            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1036            
1037    class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1038         def test_solve(self):
1039            domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
1040            x=Solution(domain).getX()
1041            # --- set exact solution ----
1042            u_ex=Vector(0,Solution(domain))
1043            u_ex[0]=1.+2.*x[0]+3.*x[1]
1044            u_ex[1]=-1.+3.*x[0]+2.*x[1]
1045            # --- set exact gradient -----------
1046            g_ex=Data(0.,(2,2),Solution(domain))
1047            g_ex[0,0]=2.
1048            g_ex[0,1]=3.
1049            g_ex[1,0]=3.
1050            g_ex[1,1]=2.
1051            # -------- test gradient --------------------------------
1052            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1053            # -------- set-up PDE -----------------------------------
1054            pde=LinearPDE(domain,numEquations=2)
1055            mask=whereZero(x[0])
1056            pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1057            A=Tensor4(0,Function(domain))
1058            A[0,:,0,:]=kronecker(2)
1059            A[1,:,1,:]=kronecker(2)
1060            Y=Vector(0.,Function(domain))
1061            Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1062            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1063            pde.setValue(A=A,
1064                         D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1065                         Y=Y,
1066                         y=matrixmult(g_ex,domain.getNormal()))
1067            # -------- get the solution ---------------------------
1068            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1069            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1070            pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1071            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1072            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1073            u=pde.getSolution()
1074            # -------- test the solution ---------------------------
1075            error=Lsup(u-u_ex)
1076            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1077    
1078    class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1079         def test_solve(self):
1080            domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
1081            x=Solution(domain).getX()
1082            # --- set exact solution ----
1083            u_ex=Vector(0,Solution(domain))
1084            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1085            u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
1086            # --- set exact gradient -----------
1087            g_ex=Data(0.,(2,2),Solution(domain))
1088            g_ex[0,0]=2.+8.*x[0]+5.*x[1]
1089            g_ex[0,1]=3.+5.*x[0]+12.*x[1]
1090            g_ex[1,0]=4.+2.*x[0]+6.*x[1]
1091            g_ex[1,1]=2.+6.*x[0]+8.*x[1]
1092            # -------- test gradient --------------------------------
1093            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1094            # -------- set-up PDE -----------------------------------
1095            pde=LinearPDE(domain,numEquations=2)
1096            mask=whereZero(x[0])
1097            pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1098            A=Tensor4(0,Function(domain))
1099            A[0,:,0,:]=kronecker(2)
1100            A[1,:,1,:]=kronecker(2)
1101            Y=Vector(0.,Function(domain))
1102            Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1103            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1104            pde.setValue(A=A,
1105                         D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1106                         Y=Y-[20.,10.],
1107                         y=matrixmult(g_ex,domain.getNormal()))
1108            # -------- get the solution ---------------------------
1109            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1110            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1111            pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1112            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1113            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1114            u=pde.getSolution()
1115            # -------- test the solution ---------------------------
1116            error=Lsup(u-u_ex)
1117            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1118    
1119    class SimpleSolve_Brick_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1120         def test_solve(self):
1121            domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1122            x=Solution(domain).getX()
1123            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
1124            # --- set exact gradient -----------
1125            g_ex=Data(0.,(3,),Solution(domain))
1126            g_ex[0]=2.
1127            g_ex[1]=3.
1128            g_ex[2]=4.
1129            # -------- test gradient --------------------------------
1130            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1131            # -------- set-up PDE -----------------------------------
1132            pde=LinearPDE(domain,numEquations=1)
1133            mask=whereZero(x[0])
1134            pde.setValue(r=u_ex,q=mask)
1135            pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
1136            # -------- get the solution ---------------------------
1137            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1138            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1139            pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1140            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1141            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1142            u=pde.getSolution()
1143            # -------- test the solution ---------------------------
1144            error=Lsup(u-u_ex)
1145            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1146            
1147    class SimpleSolve_Brick_Order1_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1148         def test_solve(self):
1149            domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1150            x=Solution(domain).getX()
1151            # --- set exact solution ----
1152            u_ex=Vector(0,Solution(domain))
1153            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
1154            u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
1155            u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
1156            # --- set exact gradient -----------
1157            g_ex=Data(0.,(3,3),Solution(domain))
1158            g_ex[0,0]=2.
1159            g_ex[0,1]=3.
1160            g_ex[0,2]=4.
1161            g_ex[1,0]=4.
1162            g_ex[1,1]=1.
1163            g_ex[1,2]=-2.
1164            g_ex[2,0]=8.
1165            g_ex[2,1]=4.
1166            g_ex[2,2]=5.
1167            # -------- test gradient --------------------------------
1168            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1169            # -------- set-up PDE -----------------------------------
1170            pde=LinearPDE(domain,numEquations=3)
1171            mask=whereZero(x[0])
1172            pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1173            A=Tensor4(0,Function(domain))
1174            A[0,:,0,:]=kronecker(3)
1175            A[1,:,1,:]=kronecker(3)
1176            A[2,:,2,:]=kronecker(3)
1177            Y=Vector(0.,Function(domain))
1178            Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1179            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1180            Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1181            pde.setValue(A=A,
1182                         D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1183                         Y=Y,
1184                         y=matrixmult(g_ex,domain.getNormal()))
1185            # -------- get the solution ---------------------------
1186            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1187            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1188            pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1189            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1190            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1191            u=pde.getSolution()
1192            # -------- test the solution ---------------------------
1193            error=Lsup(u-u_ex)
1194            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1195            
1196    class SimpleSolve_Brick_Order2_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1197         def test_solve(self):
1198            domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1199            x=Solution(domain).getX()
1200            # --- set exact solution ----
1201            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
1202            # --- set exact gradient -----------
1203            g_ex=Data(0.,(3,),Solution(domain))
1204            g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1205            g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1206            g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1207            # -------- test gradient --------------------------------
1208            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1209            # -------- set-up PDE -----------------------------------
1210            pde=LinearPDE(domain,numEquations=1)
1211            mask=whereZero(x[0])
1212            pde.setValue(r=u_ex,q=mask)
1213            pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
1214            # -------- get the solution ---------------------------
1215            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1216            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1217            pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1218            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1219            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1220            u=pde.getSolution()
1221            # -------- test the solution ---------------------------
1222            error=Lsup(u-u_ex)
1223            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1224            
1225    class SimpleSolve_Brick_Order2_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1226         def test_solve(self):
1227            domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1228            x=Solution(domain).getX()
1229            # --- set exact solution ----
1230            u_ex=Vector(0,Solution(domain))
1231            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
1232            u_ex[1]=2.+4.*x[0]+1.*x[1]-6.*x[2]+3.*x[0]*x[1]+2.*x[1]*x[2]-8.*x[2]*x[0]-2.*x[0]**2+7.*x[1]**2+5.*x[2]**2
1233            u_ex[2]=-2.+7.*x[0]+9.*x[1]+2*x[2]-6.*x[0]*x[1]+8.*x[1]*x[2]+2.*x[2]*x[0]+2.*x[0]**2+8.*x[1]**2+1.*x[2]**2
1234            # --- set exact gradient -----------
1235            g_ex=Data(0.,(3,3),Solution(domain))
1236            g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1237            g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1238            g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1239            g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
1240            g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
1241            g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
1242            g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
1243            g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
1244            g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
1245            # -------- test gradient --------------------------------
1246            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1247            # -------- set-up PDE -----------------------------------
1248            pde=LinearPDE(domain,numEquations=3)
1249            mask=whereZero(x[0])
1250            pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1251            Y=Vector(0.,Function(domain))
1252            Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1253            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1254            Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1255            A=Tensor4(0,Function(domain))
1256            A[0,:,0,:]=kronecker(3)
1257            A[1,:,1,:]=kronecker(3)
1258            A[2,:,2,:]=kronecker(3)
1259            pde.setValue(A=A,
1260                         D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1261                         Y=Y-numpy.array([60.,20.,22.]),
1262                         y=matrixmult(g_ex,domain.getNormal()))
1263            # -------- get the solution ---------------------------
1264            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1265            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1266            pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1267            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1268            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1269            u=pde.getSolution()
1270            # -------- test the solution ---------------------------
1271            error=Lsup(u-u_ex)
1272            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1273    
1274    
1275    class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
1276         def test_solve(self):
1277        # Tell about how many MPI CPUs and OpenMP threads
1278            domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
1279            x=Solution(domain).getX()
1280            # --- set exact solution ----
1281            u_ex=Scalar(0,Solution(domain))
1282            u_ex=1.+2.*x[0]+3.*x[1]
1283            # --- set exact gradient -----------
1284            g_ex=Data(0.,(2,),Solution(domain))
1285            g_ex[0]=2.
1286            g_ex[1]=3.
1287            # -------- test gradient --------------------------------
1288            g=grad(u_ex)
1289            self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
1290            # -------- set-up PDE -----------------------------------
1291            pde=LinearPDE(domain,numEquations=1)
1292            mask=whereZero(x[0])
1293            pde.setValue(r=u_ex,q=mask)
1294            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
1295            # -------- get the solution ---------------------------
1296            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1297            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1298            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1299            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1300            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1301            u=pde.getSolution()
1302            # -------- test the solution ---------------------------
1303            error=Lsup(u-u_ex)
1304            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1305    class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1306         def test_solve(self):
1307            domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
1308            x=Solution(domain).getX()
1309            # --- set exact solution ----
1310            u_ex=Scalar(0,Solution(domain))
1311            u_ex=1.+2.*x[0]+3.*x[1]
1312            # --- set exact gradient -----------
1313            g_ex=Data(0.,(2,),Solution(domain))
1314            g_ex[0]=2.
1315            g_ex[1]=3.
1316            # -------- test gradient --------------------------------
1317            g=grad(u_ex)
1318            self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
1319            # -------- set-up PDE -----------------------------------
1320            pde=LinearPDE(domain,numEquations=1)
1321            mask=whereZero(x[0])
1322            pde.setValue(r=u_ex,q=mask)
1323            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
1324            # -------- get the solution ---------------------------
1325            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1326            pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1327            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1328            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1329            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1330            u=pde.getSolution()
1331            # -------- test the solution ---------------------------
1332            error=Lsup(u-u_ex)
1333            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1334    class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1335         def test_solve(self):
1336            domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
1337            x=Solution(domain).getX()
1338            # --- set exact solution ----
1339            u_ex=Vector(0,Solution(domain))
1340            u_ex[0]=1.+2.*x[0]+3.*x[1]
1341            u_ex[1]=-1.+3.*x[0]+2.*x[1]
1342            # --- set exact gradient -----------
1343            g_ex=Data(0.,(2,2),Solution(domain))
1344            g_ex[0,0]=2.
1345            g_ex[0,1]=3.
1346            g_ex[1,0]=3.
1347            g_ex[1,1]=2.
1348            # -------- test gradient --------------------------------
1349            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1350            # -------- set-up PDE -----------------------------------
1351            pde=LinearPDE(domain,numEquations=2)
1352            mask=whereZero(x[0])
1353            pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1354            A=Tensor4(0,Function(domain))
1355            A[0,:,0,:]=kronecker(2)
1356            A[1,:,1,:]=kronecker(2)
1357            Y=Vector(0.,Function(domain))
1358            Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1359            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1360            pde.setValue(A=A,
1361                         D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1362                         Y=Y,
1363                         y=matrixmult(g_ex,domain.getNormal()))
1364            # -------- get the solution ---------------------------
1365            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1366            pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1367            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1368            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1369            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1370            u=pde.getSolution()
1371            # -------- test the solution ---------------------------
1372            error=Lsup(u-u_ex)
1373            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1374    class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1375         def test_solve(self):
1376            domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
1377            x=Solution(domain).getX()
1378            # --- set exact solution ----
1379            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1380            # --- set exact gradient -----------
1381            g_ex=Data(0.,(2,),Solution(domain))
1382            g_ex[0]=2.+8.*x[0]+5.*x[1]
1383            g_ex[1]=3.+5.*x[0]+12.*x[1]
1384            # -------- test gradient --------------------------------
1385            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1386            # -------- set-up PDE -----------------------------------
1387            pde=LinearPDE(domain,numEquations=1)
1388            mask=whereZero(x[0])
1389            pde.setValue(r=u_ex,q=mask)
1390            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
1391            # -------- get the solution ---------------------------
1392            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1393            pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1394            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1395            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1396            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1397            u=pde.getSolution()
1398            # -------- test the solution ---------------------------
1399            error=Lsup(u-u_ex)
1400            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1401    class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1402         def test_solve(self):
1403            domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
1404            x=Solution(domain).getX()
1405            # --- set exact solution ----
1406            u_ex=Vector(0,Solution(domain))
1407            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1408            u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
1409            # --- set exact gradient -----------
1410            g_ex=Data(0.,(2,2),Solution(domain))
1411            g_ex[0,0]=2.+8.*x[0]+5.*x[1]
1412            g_ex[0,1]=3.+5.*x[0]+12.*x[1]
1413            g_ex[1,0]=4.+2.*x[0]+6.*x[1]
1414            g_ex[1,1]=2.+6.*x[0]+8.*x[1]
1415            # -------- test gradient --------------------------------
1416            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1417            # -------- set-up PDE -----------------------------------
1418            pde=LinearPDE(domain,numEquations=2)
1419            mask=whereZero(x[0])
1420            pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1421            A=Tensor4(0,Function(domain))
1422            A[0,:,0,:]=kronecker(2)
1423            A[1,:,1,:]=kronecker(2)
1424            Y=Vector(0.,Function(domain))
1425            Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1426            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1427            pde.setValue(A=A,
1428                         D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1429                         Y=Y-[20.,10.],
1430                         y=matrixmult(g_ex,domain.getNormal()))
1431            # -------- get the solution ---------------------------
1432            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1433            pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1434            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1435            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1436            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1437            u=pde.getSolution()
1438            # -------- test the solution ---------------------------
1439            error=Lsup(u-u_ex)
1440            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1441    class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1442         def test_solve(self):
1443            domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1444            x=Solution(domain).getX()
1445            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
1446            # --- set exact gradient -----------
1447            g_ex=Data(0.,(3,),Solution(domain))
1448            g_ex[0]=2.
1449            g_ex[1]=3.
1450            g_ex[2]=4.
1451            # -------- test gradient --------------------------------
1452            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1453            # -------- set-up PDE -----------------------------------
1454            pde=LinearPDE(domain,numEquations=1)
1455            mask=whereZero(x[0])
1456            pde.setValue(r=u_ex,q=mask)
1457            pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
1458            # -------- get the solution ---------------------------
1459            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1460            pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1461            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1462            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1463            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1464            u=pde.getSolution()
1465            # -------- test the solution ---------------------------
1466            error=Lsup(u-u_ex)
1467            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1468    class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1469         def test_solve(self):
1470            domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1471            x=Solution(domain).getX()
1472            # --- set exact solution ----
1473            u_ex=Vector(0,Solution(domain))
1474            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
1475            u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
1476            u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
1477            # --- set exact gradient -----------
1478            g_ex=Data(0.,(3,3),Solution(domain))
1479            g_ex[0,0]=2.
1480            g_ex[0,1]=3.
1481            g_ex[0,2]=4.
1482            g_ex[1,0]=4.
1483            g_ex[1,1]=1.
1484            g_ex[1,2]=-2.
1485            g_ex[2,0]=8.
1486            g_ex[2,1]=4.
1487            g_ex[2,2]=5.
1488            # -------- test gradient --------------------------------
1489            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1490            # -------- set-up PDE -----------------------------------
1491            pde=LinearPDE(domain,numEquations=3)
1492            mask=whereZero(x[0])
1493            pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1494            A=Tensor4(0,Function(domain))
1495            A[0,:,0,:]=kronecker(3)
1496            A[1,:,1,:]=kronecker(3)
1497            A[2,:,2,:]=kronecker(3)
1498            Y=Vector(0.,Function(domain))
1499            Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1500            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1501            Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1502            pde.setValue(A=A,
1503                         D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1504                         Y=Y,
1505                         y=matrixmult(g_ex,domain.getNormal()))
1506            # -------- get the solution ---------------------------
1507            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1508            pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1509            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1510            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1511            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1512            u=pde.getSolution()
1513            # -------- test the solution ---------------------------
1514            error=Lsup(u-u_ex)
1515            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1516    class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1517         def test_solve(self):
1518            domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1519            x=Solution(domain).getX()
1520            # --- set exact solution ----
1521            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
1522            # --- set exact gradient -----------
1523            g_ex=Data(0.,(3,),Solution(domain))
1524            g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1525            g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1526            g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1527            # -------- test gradient --------------------------------
1528            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1529            # -------- set-up PDE -----------------------------------
1530            pde=LinearPDE(domain,numEquations=1)
1531            mask=whereZero(x[0])
1532            pde.setValue(r=u_ex,q=mask)
1533            pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
1534            # -------- get the solution ---------------------------
1535            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1536            pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1537            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1538            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1539            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1540            u=pde.getSolution()
1541            # -------- test the solution ---------------------------
1542            error=Lsup(u-u_ex)
1543            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1544    class SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1545         def test_solve(self):
1546            domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1547            x=Solution(domain).getX()
1548            # --- set exact solution ----
1549            u_ex=Vector(0,Solution(domain))
1550            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
1551            u_ex[1]=2.+4.*x[0]+1.*x[1]-6.*x[2]+3.*x[0]*x[1]+2.*x[1]*x[2]-8.*x[2]*x[0]-2.*x[0]**2+7.*x[1]**2+5.*x[2]**2
1552            u_ex[2]=-2.+7.*x[0]+9.*x[1]+2*x[2]-6.*x[0]*x[1]+8.*x[1]*x[2]+2.*x[2]*x[0]+2.*x[0]**2+8.*x[1]**2+1.*x[2]**2
1553            # --- set exact gradient -----------
1554            g_ex=Data(0.,(3,3),Solution(domain))
1555            g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1556            g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1557            g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1558            g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
1559            g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
1560            g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
1561            g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
1562            g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
1563            g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
1564            # -------- test gradient --------------------------------
1565            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1566            # -------- set-up PDE -----------------------------------
1567            pde=LinearPDE(domain,numEquations=3)
1568            mask=whereZero(x[0])
1569            pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1570            Y=Vector(0.,Function(domain))
1571            Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1572            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1573            Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1574            A=Tensor4(0,Function(domain))
1575            A[0,:,0,:]=kronecker(3)
1576            A[1,:,1,:]=kronecker(3)
1577            A[2,:,2,:]=kronecker(3)
1578            pde.setValue(A=A,
1579                         D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1580                         Y=Y-numpy.array([60.,20.,22.]),
1581                         y=matrixmult(g_ex,domain.getNormal()))
1582            # -------- get the solution ---------------------------
1583            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1584            pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1585            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1586            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1587            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1588            u=pde.getSolution()
1589            # -------- test the solution ---------------------------
1590            error=Lsup(u-u_ex)
1591            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1592    
1593    class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1594         def test_solve(self):
1595            domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
1596            x=Solution(domain).getX()
1597            # --- set exact solution ----
1598            u_ex=Scalar(0,Solution(domain))
1599            u_ex=1.+2.*x[0]+3.*x[1]
1600            # --- set exact gradient -----------
1601            g_ex=Data(0.,(2,),Solution(domain))
1602            g_ex[0]=2.
1603            g_ex[1]=3.
1604            # -------- test gradient --------------------------------
1605            g=grad(u_ex)
1606            self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
1607            # -------- set-up PDE -----------------------------------
1608            pde=LinearPDE(domain,numEquations=1)
1609            mask=whereZero(x[0])
1610            pde.setValue(r=u_ex,q=mask)
1611            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
1612            # -------- get the solution ---------------------------
1613            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1614            pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1615            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1616            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1617            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1618            u=pde.getSolution()
1619            # -------- test the solution ---------------------------
1620            error=Lsup(u-u_ex)
1621            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1622    
1623    class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1624         def test_solve(self):
1625            domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
1626            x=Solution(domain).getX()
1627            # --- set exact solution ----
1628            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1629            # --- set exact gradient -----------
1630            g_ex=Data(0.,(2,),Solution(domain))
1631            g_ex[0]=2.+8.*x[0]+5.*x[1]
1632            g_ex[1]=3.+5.*x[0]+12.*x[1]
1633            # -------- test gradient --------------------------------
1634            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1635            # -------- set-up PDE -----------------------------------
1636            pde=LinearPDE(domain,numEquations=1)
1637            mask=whereZero(x[0])
1638            pde.setValue(r=u_ex,q=mask)
1639            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
1640            # -------- get the solution ---------------------------
1641            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1642            pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1643            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1644            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1645            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1646            u=pde.getSolution()
1647            # -------- test the solution ---------------------------
1648            error=Lsup(u-u_ex)
1649            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1650            
1651    class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1652         def test_solve(self):
1653            domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
1654            x=Solution(domain).getX()
1655            # --- set exact solution ----
1656            u_ex=Vector(0,Solution(domain))
1657            u_ex[0]=1.+2.*x[0]+3.*x[1]
1658            u_ex[1]=-1.+3.*x[0]+2.*x[1]
1659            # --- set exact gradient -----------
1660            g_ex=Data(0.,(2,2),Solution(domain))
1661            g_ex[0,0]=2.
1662            g_ex[0,1]=3.
1663            g_ex[1,0]=3.
1664            g_ex[1,1]=2.
1665            # -------- test gradient --------------------------------
1666            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1667            # -------- set-up PDE -----------------------------------
1668            pde=LinearPDE(domain,numEquations=2)
1669            mask=whereZero(x[0])
1670            pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1671            A=Tensor4(0,Function(domain))
1672            A[0,:,0,:]=kronecker(2)
1673            A[1,:,1,:]=kronecker(2)
1674            Y=Vector(0.,Function(domain))
1675            Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1676            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1677            pde.setValue(A=A,
1678                         D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1679                         Y=Y,
1680                         y=matrixmult(g_ex,domain.getNormal()))
1681            # -------- get the solution ---------------------------
1682            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1683            pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1684            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1685            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1686            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1687            u=pde.getSolution()
1688            # -------- test the solution ---------------------------
1689            error=Lsup(u-u_ex)
1690            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1691    
1692    class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1693         def test_solve(self):
1694            domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
1695            x=Solution(domain).getX()
1696            # --- set exact solution ----
1697            u_ex=Vector(0,Solution(domain))
1698            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1699            u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
1700            # --- set exact gradient -----------
1701            g_ex=Data(0.,(2,2),Solution(domain))
1702            g_ex[0,0]=2.+8.*x[0]+5.*x[1]
1703            g_ex[0,1]=3.+5.*x[0]+12.*x[1]
1704            g_ex[1,0]=4.+2.*x[0]+6.*x[1]
1705            g_ex[1,1]=2.+6.*x[0]+8.*x[1]
1706            # -------- test gradient --------------------------------
1707            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1708            # -------- set-up PDE -----------------------------------
1709            pde=LinearPDE(domain,numEquations=2)
1710            mask=whereZero(x[0])
1711            pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1712            A=Tensor4(0,Function(domain))
1713            A[0,:,0,:]=kronecker(2)
1714            A[1,:,1,:]=kronecker(2)
1715            Y=Vector(0.,Function(domain))
1716            Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1717            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1718            pde.setValue(A=A,
1719                         D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1720                         Y=Y-[20.,10.],
1721                         y=matrixmult(g_ex,domain.getNormal()))
1722            # -------- get the solution ---------------------------
1723            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1724            pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1725            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1726            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1727            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1728            u=pde.getSolution()
1729            # -------- test the solution ---------------------------
1730            error=Lsup(u-u_ex)
1731            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1732    
1733    class SimpleSolve_Brick_Order1_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1734         def test_solve(self):
1735            domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1736            x=Solution(domain).getX()
1737            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
1738            # --- set exact gradient -----------
1739            g_ex=Data(0.,(3,),Solution(domain))
1740            g_ex[0]=2.
1741            g_ex[1]=3.
1742            g_ex[2]=4.
1743            # -------- test gradient --------------------------------
1744            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1745            # -------- set-up PDE -----------------------------------
1746            pde=LinearPDE(domain,numEquations=1)
1747            mask=whereZero(x[0])
1748            pde.setValue(r=u_ex,q=mask)
1749            pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
1750            # -------- get the solution ---------------------------
1751            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1752            pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1753            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1754            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1755            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1756            u=pde.getSolution()
1757            # -------- test the solution ---------------------------
1758            error=Lsup(u-u_ex)
1759            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1760            
1761    class SimpleSolve_Brick_Order1_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1762         def test_solve(self):
1763            domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1764            x=Solution(domain).getX()
1765            # --- set exact solution ----
1766            u_ex=Vector(0,Solution(domain))
1767            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
1768            u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
1769            u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
1770            # --- set exact gradient -----------
1771            g_ex=Data(0.,(3,3),Solution(domain))
1772            g_ex[0,0]=2.
1773            g_ex[0,1]=3.
1774            g_ex[0,2]=4.
1775            g_ex[1,0]=4.
1776            g_ex[1,1]=1.
1777            g_ex[1,2]=-2.
1778            g_ex[2,0]=8.
1779            g_ex[2,1]=4.
1780            g_ex[2,2]=5.
1781            # -------- test gradient --------------------------------
1782            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1783            # -------- set-up PDE -----------------------------------
1784            pde=LinearPDE(domain,numEquations=3)
1785            mask=whereZero(x[0])
1786            pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1787            A=Tensor4(0,Function(domain))
1788            A[0,:,0,:]=kronecker(3)
1789            A[1,:,1,:]=kronecker(3)
1790            A[2,:,2,:]=kronecker(3)
1791            Y=Vector(0.,Function(domain))
1792            Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1793            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1794            Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1795            pde.setValue(A=A,
1796                         D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1797                         Y=Y,
1798                         y=matrixmult(g_ex,domain.getNormal()))
1799            # -------- get the solution ---------------------------
1800            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1801            pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1802            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1803            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1804            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1805            u=pde.getSolution()
1806            # -------- test the solution ---------------------------
1807            error=Lsup(u-u_ex)
1808            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1809            
1810    class SimpleSolve_Brick_Order2_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1811         def test_solve(self):
1812            domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1813            x=Solution(domain).getX()
1814            # --- set exact solution ----
1815            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
1816            # --- set exact gradient -----------
1817            g_ex=Data(0.,(3,),Solution(domain))
1818            g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1819            g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1820            g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1821            # -------- test gradient --------------------------------
1822            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1823            # -------- set-up PDE -----------------------------------
1824            pde=LinearPDE(domain,numEquations=1)
1825            mask=whereZero(x[0])
1826            pde.setValue(r=u_ex,q=mask)
1827            pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
1828            # -------- get the solution ---------------------------
1829            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1830            pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1831            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1832            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1833            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1834            u=pde.getSolution()
1835            # -------- test the solution ---------------------------
1836            error=Lsup(u-u_ex)
1837            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1838            
1839    class SimpleSolve_Brick_Order2_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1840         def test_solve(self):
1841            domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1842            x=Solution(domain).getX()
1843            # --- set exact solution ----
1844            u_ex=Vector(0,Solution(domain))
1845            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
1846            u_ex[1]=2.+4.*x[0]+1.*x[1]-6.*x[2]+3.*x[0]*x[1]+2.*x[1]*x[2]-8.*x[2]*x[0]-2.*x[0]**2+7.*x[1]**2+5.*x[2]**2
1847            u_ex[2]=-2.+7.*x[0]+9.*x[1]+2*x[2]-6.*x[0]*x[1]+8.*x[1]*x[2]+2.*x[2]*x[0]+2.*x[0]**2+8.*x[1]**2+1.*x[2]**2
1848            # --- set exact gradient -----------
1849            g_ex=Data(0.,(3,3),Solution(domain))
1850            g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1851            g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1852            g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1853            g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
1854            g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
1855            g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
1856            g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
1857            g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
1858            g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
1859            # -------- test gradient --------------------------------
1860            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1861            # -------- set-up PDE -----------------------------------
1862            pde=LinearPDE(domain,numEquations=3)
1863            mask=whereZero(x[0])
1864            pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1865            Y=Vector(0.,Function(domain))
1866            Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1867            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1868            Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1869            A=Tensor4(0,Function(domain))
1870            A[0,:,0,:]=kronecker(3)
1871            A[1,:,1,:]=kronecker(3)
1872            A[2,:,2,:]=kronecker(3)
1873            pde.setValue(A=A,
1874                         D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1875                         Y=Y-numpy.array([60.,20.,22.]),
1876                         y=matrixmult(g_ex,domain.getNormal()))
1877            # -------- get the solution ---------------------------
1878            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1879            pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1880            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1881            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1882            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1883            u=pde.getSolution()
1884            # -------- test the solution ---------------------------
1885            error=Lsup(u-u_ex)
1886            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1887    
1888    
1889    class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_MINRES_GaussSeidel(unittest.TestCase):
1890         def test_solve(self):
1891            domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
1892            x=Solution(domain).getX()
1893            # --- set exact solution ----
1894            u_ex=Scalar(0,Solution(domain))
1895            u_ex=1.+2.*x[0]+3.*x[1]
1896            # --- set exact gradient -----------
1897            g_ex=Data(0.,(2,),Solution(domain))
1898            g_ex[0]=2.
1899            g_ex[1]=3.
1900            # -------- test gradient --------------------------------
1901            g=grad(u_ex)
1902            self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
1903            # -------- set-up PDE -----------------------------------
1904            pde=LinearPDE(domain,numEquations=1)
1905            mask=whereZero(x[0])
1906            pde.setValue(r=u_ex,q=mask)
1907            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
1908            # -------- get the solution ---------------------------
1909            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1910            pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
1911            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1912            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1913            
1914            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1915            u=pde.getSolution()
1916            # -------- test the solution ---------------------------
1917            error=Lsup(u-u_ex)
1918            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1919    
1920    class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_MINRES_GaussSeidel(unittest.TestCase):
1921         def test_solve(self):
1922            domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
1923            x=Solution(domain).getX()
1924            # --- set exact solution ----
1925            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1926            # --- set exact gradient -----------
1927            g_ex=Data(0.,(2,),Solution(domain))
1928            g_ex[0]=2.+8.*x[0]+5.*x[1]
1929            g_ex[1]=3.+5.*x[0]+12.*x[1]
1930            # -------- test gradient --------------------------------
1931            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1932            # -------- set-up PDE -----------------------------------
1933            pde=LinearPDE(domain,numEquations=1)
1934            mask=whereZero(x[0])
1935            pde.setValue(r=u_ex,q=mask)
1936            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
1937            # -------- get the solution ---------------------------
1938            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1939            pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
1940            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1941            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1942            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1943            u=pde.getSolution()
1944            # -------- test the solution ---------------------------
1945            error=Lsup(u-u_ex)
1946            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1947            
1948    class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_MINRES_GaussSeidel(unittest.TestCase):
1949         def test_solve(self):
1950            domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
1951            x=Solution(domain).getX()
1952            # --- set exact solution ----
1953            u_ex=Vector(0,Solution(domain))
1954            u_ex[0]=1.+2.*x[0]+3.*x[1]
1955            u_ex[1]=-1.+3.*x[0]+2.*x[1]
1956            # --- set exact gradient -----------
1957            g_ex=Data(0.,(2,2),Solution(domain))
1958            g_ex[0,0]=2.
1959            g_ex[0,1]=3.
1960            g_ex[1,0]=3.
1961            g_ex[1,1]=2.
1962            # -------- test gradient --------------------------------
1963            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1964            # -------- set-up PDE -----------------------------------
1965            pde=LinearPDE(domain,numEquations=2)
1966            mask=whereZero(x[0])
1967            pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1968            A=Tensor4(0,Function(domain))
1969            A[0,:,0,:]=kronecker(2)
1970            A[1,:,1,:]=kronecker(2)
1971            Y=Vector(0.,Function(domain))
1972            Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1973            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1974            pde.setValue(A=A,
1975                         D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1976                         Y=Y,
1977                         y=matrixmult(g_ex,domain.getNormal()))
1978            # -------- get the solution ---------------------------
1979            pde.getSolverOptions().setTolerance(SOLVER_TOL)
1980            pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
1981            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1982            pde.getSolverOptions().setPackage(SolverOptions.PASO)
1983            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1984            u=pde.getSolution()
1985            # -------- test the solution ---------------------------
1986            error=Lsup(u-u_ex)
1987            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1988            
1989    class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_MINRES_GaussSeidel(unittest.TestCase):
1990         def test_solve(self):
1991            domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
1992            x=Solution(domain).getX()
1993            # --- set exact solution ----
1994            u_ex=Vector(0,Solution(domain))
1995            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1996            u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
1997            # --- set exact gradient -----------
1998            g_ex=Data(0.,(2,2),Solution(domain))
1999            g_ex[0,0]=2.+8.*x[0]+5.*x[1]
2000            g_ex[0,1]=3.+5.*x[0]+12.*x[1]
2001            g_ex[1,0]=4.+2.*x[0]+6.*x[1]
2002            g_ex[1,1]=2.+6.*x[0]+8.*x[1]
2003            # -------- test gradient --------------------------------
2004            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2005            # -------- set-up PDE -----------------------------------
2006            pde=LinearPDE(domain,numEquations=2)
2007            mask=whereZero(x[0])
2008            pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
2009            A=Tensor4(0,Function(domain))
2010            A[0,:,0,:]=kronecker(2)
2011            A[1,:,1,:]=kronecker(2)
2012            Y=Vector(0.,Function(domain))
2013            Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
2014            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
2015            pde.setValue(A=A,
2016                         D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
2017                         Y=Y-[20.,10.],
2018                         y=matrixmult(g_ex,domain.getNormal()))
2019            # -------- get the solution ---------------------------
2020            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2021            pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
2022            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2023            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2024            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2025            u=pde.getSolution()
2026            # -------- test the solution ---------------------------
2027            error=Lsup(u-u_ex)
2028            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2029    
2030    class SimpleSolve_Brick_Order1_SinglePDE_Paso_MINRES_GaussSeidel(unittest.TestCase):
2031         def test_solve(self):
2032            domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
2033            x=Solution(domain).getX()
2034            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
2035            # --- set exact gradient -----------
2036            g_ex=Data(0.,(3,),Solution(domain))
2037            g_ex[0]=2.
2038            g_ex[1]=3.
2039            g_ex[2]=4.
2040            # -------- test gradient --------------------------------
2041            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2042            # -------- set-up PDE -----------------------------------
2043            pde=LinearPDE(domain,numEquations=1)
2044            mask=whereZero(x[0])
2045            pde.setValue(r=u_ex,q=mask)
2046            pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
2047            # -------- get the solution ---------------------------
2048            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2049            pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
2050            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2051            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2052            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2053            u=pde.getSolution()
2054            # -------- test the solution ---------------------------
2055            error=Lsup(u-u_ex)
2056            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2057            
2058    class SimpleSolve_Brick_Order1_SystemPDE_Paso_MINRES_GaussSeidel(unittest.TestCase):
2059         def test_solve(self):
2060            domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
2061            x=Solution(domain).getX()
2062            # --- set exact solution ----
2063            u_ex=Vector(0,Solution(domain))
2064            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
2065            u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
2066            u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
2067            # --- set exact gradient -----------
2068            g_ex=Data(0.,(3,3),Solution(domain))
2069            g_ex[0,0]=2.
2070            g_ex[0,1]=3.
2071            g_ex[0,2]=4.
2072            g_ex[1,0]=4.
2073            g_ex[1,1]=1.
2074            g_ex[1,2]=-2.
2075            g_ex[2,0]=8.
2076            g_ex[2,1]=4.
2077            g_ex[2,2]=5.
2078            # -------- test gradient --------------------------------
2079            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2080            # -------- set-up PDE -----------------------------------
2081            pde=LinearPDE(domain,numEquations=3)
2082            mask=whereZero(x[0])
2083            pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
2084            A=Tensor4(0,Function(domain))
2085            A[0,:,0,:]=kronecker(3)
2086            A[1,:,1,:]=kronecker(3)
2087            A[2,:,2,:]=kronecker(3)
2088            Y=Vector(0.,Function(domain))
2089            Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
2090            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
2091            Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
2092            pde.setValue(A=A,
2093                         D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
2094                         Y=Y,
2095                         y=matrixmult(g_ex,domain.getNormal()))
2096            # -------- get the solution ---------------------------
2097            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2098            pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
2099            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2100            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2101            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2102            u=pde.getSolution()
2103            # -------- test the solution ---------------------------
2104            error=Lsup(u-u_ex)
2105            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2106            
2107    class SimpleSolve_Brick_Order2_SinglePDE_Paso_MINRES_GaussSeidel(unittest.TestCase):
2108         def test_solve(self):
2109            domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
2110            x=Solution(domain).getX()
2111            # --- set exact solution ----
2112            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
2113            # --- set exact gradient -----------
2114            g_ex=Data(0.,(3,),Solution(domain))
2115            g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
2116            g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
2117            g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
2118            # -------- test gradient --------------------------------
2119            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2120            # -------- set-up PDE -----------------------------------
2121            pde=LinearPDE(domain,numEquations=1)
2122            mask=whereZero(x[0])
2123            pde.setValue(r=u_ex,q=mask)
2124            pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
2125            # -------- get the solution ---------------------------
2126            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2127            pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
2128            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2129            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2130            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2131            u=pde.getSolution()
2132            # -------- test the solution ---------------------------
2133            error=Lsup(u-u_ex)
2134            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2135            
2136    class SimpleSolve_Brick_Order2_SystemPDE_Paso_MINRES_GaussSeidel(unittest.TestCase):
2137         def test_solve(self):
2138            domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
2139            x=Solution(domain).getX()
2140            # --- set exact solution ----
2141            u_ex=Vector(0,Solution(domain))
2142            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
2143            u_ex[1]=2.+4.*x[0]+1.*x[1]-6.*x[2]+3.*x[0]*x[1]+2.*x[1]*x[2]-8.*x[2]*x[0]-2.*x[0]**2+7.*x[1]**2+5.*x[2]**2
2144            u_ex[2]=-2.+7.*x[0]+9.*x[1]+2*x[2]-6.*x[0]*x[1]+8.*x[1]*x[2]+2.*x[2]*x[0]+2.*x[0]**2+8.*x[1]**2+1.*x[2]**2
2145            # --- set exact gradient -----------
2146            g_ex=Data(0.,(3,3),Solution(domain))
2147            g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
2148            g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
2149            g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
2150            g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
2151            g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
2152            g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
2153            g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
2154            g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
2155            g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
2156            # -------- test gradient --------------------------------
2157            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2158            # -------- set-up PDE -----------------------------------
2159            pde=LinearPDE(domain,numEquations=3)
2160            mask=whereZero(x[0])
2161            pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
2162            Y=Vector(0.,Function(domain))
2163            Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
2164            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
2165            Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
2166            A=Tensor4(0,Function(domain))
2167            A[0,:,0,:]=kronecker(3)
2168            A[1,:,1,:]=kronecker(3)
2169            A[2,:,2,:]=kronecker(3)
2170            pde.setValue(A=A,
2171                         D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
2172                         Y=Y-numpy.array([60.,20.,22.]),
2173                         y=matrixmult(g_ex,domain.getNormal()))
2174            # -------- get the solution ---------------------------
2175            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2176            pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
2177            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2178            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2179            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2180            u=pde.getSolution()
2181            # -------- test the solution ---------------------------
2182            error=Lsup(u-u_ex)
2183            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2184    
2185    class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
2186         def test_solve(self):
2187            domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
2188            x=Solution(domain).getX()
2189            # --- set exact solution ----
2190            u_ex=Scalar(0,Solution(domain))
2191            u_ex=1.+2.*x[0]+3.*x[1]
2192            # --- set exact gradient -----------
2193            g_ex=Data(0.,(2,),Solution(domain))
2194            g_ex[0]=2.
2195            g_ex[1]=3.
2196            # -------- test gradient --------------------------------
2197            g=grad(u_ex)
2198            self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
2199            # -------- set-up PDE -----------------------------------
2200            pde=LinearPDE(domain,numEquations=1)
2201            mask=whereZero(x[0])
2202            pde.setValue(r=u_ex,q=mask)
2203            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
2204            # -------- get the solution ---------------------------
2205            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2206            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
2207            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2208            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2209            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2210            u=pde.getSolution()
2211            # -------- test the solution ---------------------------
2212            error=Lsup(u-u_ex)
2213            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2214    
2215    class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
2216         def test_solve(self):
2217            domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
2218            x=Solution(domain).getX()
2219            # --- set exact solution ----
2220            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
2221            # --- set exact gradient -----------
2222            g_ex=Data(0.,(2,),Solution(domain))
2223            g_ex[0]=2.+8.*x[0]+5.*x[1]
2224            g_ex[1]=3.+5.*x[0]+12.*x[1]
2225            # -------- test gradient --------------------------------
2226            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2227            # -------- set-up PDE -----------------------------------
2228            pde=LinearPDE(domain,numEquations=1)
2229            mask=whereZero(x[0])
2230            pde.setValue(r=u_ex,q=mask)
2231            pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
2232            # -------- get the solution ---------------------------
2233            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2234            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
2235            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2236            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2237            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2238            u=pde.getSolution()
2239            # -------- test the solution ---------------------------
2240            error=Lsup(u-u_ex)
2241            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2242            
2243    class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
2244         def test_solve(self):
2245            domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
2246            x=Solution(domain).getX()
2247            # --- set exact solution ----
2248            u_ex=Vector(0,Solution(domain))
2249            u_ex[0]=1.+2.*x[0]+3.*x[1]
2250            u_ex[1]=-1.+3.*x[0]+2.*x[1]
2251            # --- set exact gradient -----------
2252            g_ex=Data(0.,(2,2),Solution(domain))
2253            g_ex[0,0]=2.
2254            g_ex[0,1]=3.
2255            g_ex[1,0]=3.
2256            g_ex[1,1]=2.
2257            # -------- test gradient --------------------------------
2258            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2259            # -------- set-up PDE -----------------------------------
2260            pde=LinearPDE(domain,numEquations=2)
2261            mask=whereZero(x[0])
2262            pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
2263            A=Tensor4(0,Function(domain))
2264            A[0,:,0,:]=kronecker(2)
2265            A[1,:,1,:]=kronecker(2)
2266            Y=Vector(0.,Function(domain))
2267            Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
2268            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
2269            pde.setValue(A=A,
2270                         D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
2271                         Y=Y,
2272                         y=matrixmult(g_ex,domain.getNormal()))
2273            # -------- get the solution ---------------------------
2274            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2275            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
2276            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2277            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2278            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2279            u=pde.getSolution()
2280            # -------- test the solution ---------------------------
2281            error=Lsup(u-u_ex)
2282            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2283    
2284    class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
2285         def test_solve(self):
2286            domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
2287            x=Solution(domain).getX()
2288            # --- set exact solution ----
2289            u_ex=Vector(0,Solution(domain))
2290            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
2291            u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
2292            # --- set exact gradient -----------
2293            g_ex=Data(0.,(2,2),Solution(domain))
2294            g_ex[0,0]=2.+8.*x[0]+5.*x[1]
2295            g_ex[0,1]=3.+5.*x[0]+12.*x[1]
2296            g_ex[1,0]=4.+2.*x[0]+6.*x[1]
2297            g_ex[1,1]=2.+6.*x[0]+8.*x[1]
2298            # -------- test gradient --------------------------------
2299            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2300            # -------- set-up PDE -----------------------------------
2301            pde=LinearPDE(domain,numEquations=2)
2302            mask=whereZero(x[0])
2303            pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
2304            A=Tensor4(0,Function(domain))
2305            A[0,:,0,:]=kronecker(2)
2306            A[1,:,1,:]=kronecker(2)
2307            Y=Vector(0.,Function(domain))
2308            Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
2309            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
2310            pde.setValue(A=A,
2311                         D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
2312                         Y=Y-[20.,10.],
2313                         y=matrixmult(g_ex,domain.getNormal()))
2314            # -------- get the solution ---------------------------
2315            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2316            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
2317            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2318            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2319            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2320            u=pde.getSolution()
2321            # -------- test the solution ---------------------------
2322            error=Lsup(u-u_ex)
2323            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2324    
2325    class SimpleSolve_Brick_Order1_SinglePDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
2326         def test_solve(self):
2327            domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
2328            x=Solution(domain).getX()
2329            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
2330            # --- set exact gradient -----------
2331            g_ex=Data(0.,(3,),Solution(domain))
2332            g_ex[0]=2.
2333            g_ex[1]=3.
2334            g_ex[2]=4.
2335            # -------- test gradient --------------------------------
2336            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2337            # -------- set-up PDE -----------------------------------
2338            pde=LinearPDE(domain,numEquations=1)
2339            mask=whereZero(x[0])
2340            pde.setValue(r=u_ex,q=mask)
2341            pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
2342            # -------- get the solution ---------------------------
2343            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2344            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
2345            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2346            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2347            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2348            u=pde.getSolution()
2349            # -------- test the solution ---------------------------
2350            error=Lsup(u-u_ex)
2351            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2352            
2353    class SimpleSolve_Brick_Order1_SystemPDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
2354         def test_solve(self):
2355            domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
2356            x=Solution(domain).getX()
2357            # --- set exact solution ----
2358            u_ex=Vector(0,Solution(domain))
2359            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
2360            u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
2361            u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
2362            # --- set exact gradient -----------
2363            g_ex=Data(0.,(3,3),Solution(domain))
2364            g_ex[0,0]=2.
2365            g_ex[0,1]=3.
2366            g_ex[0,2]=4.
2367            g_ex[1,0]=4.
2368            g_ex[1,1]=1.
2369            g_ex[1,2]=-2.
2370            g_ex[2,0]=8.
2371            g_ex[2,1]=4.
2372            g_ex[2,2]=5.
2373            # -------- test gradient --------------------------------
2374            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2375            # -------- set-up PDE -----------------------------------
2376            pde=LinearPDE(domain,numEquations=3)
2377            mask=whereZero(x[0])
2378            pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
2379            A=Tensor4(0,Function(domain))
2380            A[0,:,0,:]=kronecker(3)
2381            A[1,:,1,:]=kronecker(3)
2382            A[2,:,2,:]=kronecker(3)
2383            Y=Vector(0.,Function(domain))
2384            Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
2385            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
2386            Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
2387            pde.setValue(A=A,
2388                         D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
2389                         Y=Y,
2390                         y=matrixmult(g_ex,domain.getNormal()))
2391            # -------- get the solution ---------------------------
2392            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2393            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
2394            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2395            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2396            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2397            u=pde.getSolution()
2398            # -------- test the solution ---------------------------
2399            error=Lsup(u-u_ex)
2400            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2401            
2402    class SimpleSolve_Brick_Order2_SinglePDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
2403         def test_solve(self):
2404            domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
2405            x=Solution(domain).getX()
2406            # --- set exact solution ----
2407            u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
2408            # --- set exact gradient -----------
2409            g_ex=Data(0.,(3,),Solution(domain))
2410            g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
2411            g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
2412            g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
2413            # -------- test gradient --------------------------------
2414            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2415            # -------- set-up PDE -----------------------------------
2416            pde=LinearPDE(domain,numEquations=1)
2417            mask=whereZero(x[0])
2418            pde.setValue(r=u_ex,q=mask)
2419            pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
2420            # -------- get the solution ---------------------------
2421            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2422            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
2423            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2424            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2425            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2426            u=pde.getSolution()
2427            # -------- test the solution ---------------------------
2428            error=Lsup(u-u_ex)
2429            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2430            
2431    class SimpleSolve_Brick_Order2_SystemPDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
2432         def test_solve(self):
2433            domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
2434            x=Solution(domain).getX()
2435            # --- set exact solution ----
2436            u_ex=Vector(0,Solution(domain))
2437            u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
2438            u_ex[1]=2.+4.*x[0]+1.*x[1]-6.*x[2]+3.*x[0]*x[1]+2.*x[1]*x[2]-8.*x[2]*x[0]-2.*x[0]**2+7.*x[1]**2+5.*x[2]**2
2439            u_ex[2]=-2.+7.*x[0]+9.*x[1]+2*x[2]-6.*x[0]*x[1]+8.*x[1]*x[2]+2.*x[2]*x[0]+2.*x[0]**2+8.*x[1]**2+1.*x[2]**2
2440            # --- set exact gradient -----------
2441            g_ex=Data(0.,(3,3),Solution(domain))
2442            g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
2443            g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
2444            g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
2445            g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
2446            g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
2447            g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
2448            g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
2449            g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
2450            g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
2451            # -------- test gradient --------------------------------
2452            self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
2453            # -------- set-up PDE -----------------------------------
2454            pde=LinearPDE(domain,numEquations=3)
2455            mask=whereZero(x[0])
2456            pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
2457            Y=Vector(0.,Function(domain))
2458            Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
2459            Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
2460            Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
2461            A=Tensor4(0,Function(domain))
2462            A[0,:,0,:]=kronecker(3)
2463            A[1,:,1,:]=kronecker(3)
2464            A[2,:,2,:]=kronecker(3)
2465            pde.setValue(A=A,
2466                         D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
2467                         Y=Y-numpy.array([60.,20.,22.]),
2468                         y=matrixmult(g_ex,domain.getNormal()))
2469            # -------- get the solution ---------------------------
2470            pde.getSolverOptions().setTolerance(SOLVER_TOL)
2471            pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
2472            pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
2473            pde.getSolverOptions().setPackage(SolverOptions.PASO)
2474            pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
2475            u=pde.getSolution()
2476            # -------- test the solution ---------------------------
2477            error=Lsup(u-u_ex)
2478            self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
2479    
2480    
2481  if __name__ == '__main__':  if __name__ == '__main__':
2482     suite = unittest.TestSuite()     suite = unittest.TestSuite()
    suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi))  
2483     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi))     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi))
2484     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi))     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi))
2485     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi))     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi))
# Line 1003  if __name__ == '__main__': Line 2504  if __name__ == '__main__':
2504     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_MINRES_Jacobi))     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_MINRES_Jacobi))
2505     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_MINRES_Jacobi))     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_MINRES_Jacobi))
2506     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_MINRES_Jacobi))     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_MINRES_Jacobi))
2507       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi))
2508       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_BICGSTAB_Jacobi))
2509       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi))
2510       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_BICGSTAB_Jacobi))
2511       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_BICGSTAB_Jacobi))
2512       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_BICGSTAB_Jacobi))
2513       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_BICGSTAB_Jacobi))
2514       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_BICGSTAB_Jacobi))
2515    
2516       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_GaussSeidel))
2517       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_GaussSeidel))
2518       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_GaussSeidel))
2519       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_GaussSeidel))
2520       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_GaussSeidel))
2521       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_GaussSeidel))
2522       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_GaussSeidel))
2523       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_GaussSeidel))
2524       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_TFQMR_GaussSeidel))
2525       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_TFQMR_GaussSeidel))
2526       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_TFQMR_GaussSeidel))
2527       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_TFQMR_GaussSeidel))
2528       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_TFQMR_GaussSeidel))
2529       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_TFQMR_GaussSeidel))
2530       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_TFQMR_GaussSeidel))
2531       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_TFQMR_GaussSeidel))
2532       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_MINRES_GaussSeidel))
2533       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_MINRES_GaussSeidel))
2534       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_MINRES_GaussSeidel))
2535       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_MINRES_GaussSeidel))
2536       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_MINRES_GaussSeidel))
2537       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_MINRES_GaussSeidel))
2538       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_MINRES_GaussSeidel))
2539       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_MINRES_GaussSeidel))
2540       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_GaussSeidel))
2541       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_BICGSTAB_GaussSeidel))
2542       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_GaussSeidel))
2543       suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_BICGSTAB_GaussSeidel))
2544       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_BICGSTAB_GaussSeidel))
2545       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_BICGSTAB_GaussSeidel))
2546       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_BICGSTAB_GaussSeidel))
2547       suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_BICGSTAB_GaussSeidel))
2548    
2549     s=unittest.TextTestRunner(verbosity=2).run(suite)     s=unittest.TextTestRunner(verbosity=2).run(suite)
2550     if not s.wasSuccessful(): sys.exit(1)     if not s.wasSuccessful(): sys.exit(1)

Legend:
Removed from v.3096  
changed lines
  Added in v.3097

  ViewVC Help
Powered by ViewVC 1.1.26