debug

Author: cgcgcg

base/PyNucleus_base/solvers.pyx

@@ -154,6 +154,7 @@ cdef class lu_solver(solver):
             perm_r = self.perm_r
             perm_c = self.perm_c
             try:
+                print('lu1 pre', np.linalg.norm(np.asarray(b)))
                 temp = self.temp_mem
                 n = perm_c.shape[0]
                 for i in range(n):
@@ -166,14 +167,19 @@ cdef class lu_solver(solver):
                 backward_solve_csc(self.U.indptr, self.U.indices, self.U.data, x, temp)
                 for i in range(n):
                     x[i] = temp[perm_c[i]]
+                print('lu1 post', np.linalg.norm(np.asarray(x)))
             except AttributeError:
+                print('lu2 pre', np.linalg.norm(np.asarray(b)))
                 x[:] = self.Ainv.solve(np.array(b, copy=False, dtype=REAL))
+                print('lu2 post', np.linalg.norm(np.asarray(x)))
         else:
             from scipy.linalg import lu_solve
+            print('lu3 pre', np.linalg.norm(np.asarray(b)))
             assign(x, b)
             lu_solve((self.lu, self.perm),
                      np.array(x, copy=False, dtype=REAL),
                      overwrite_b=True)
+            print('lu3 post', np.linalg.norm(np.asarray(x)))
         return 1
 
     def __str__(self):

drivers/runParallelGMG.py

@@ -166,6 +166,7 @@
 for cycle, label in [(V, 'MG'),
                      (FMG_V, 'FMG')]:
     if getattr(d, 'do'+label):
+        assert np.all(np.isfinite(rhs.toarray())), rhs.toarray()
         ml.cycle = cycle
         with d.timer('Solve '+label):
             numIter = ml(rhs, x)
@@ -176,6 +177,7 @@
         its.add('Number of iterations '+label, numIter)
         res.add('Residual norm '+label, resNorm)
         resHist.add(label, residuals)
+        print(label, r0, residuals, numIter, resNorm)
 
 # set up cg
 cg = solverFactory.build('cg', A=A, maxIter=d.maxiter, tolerance=tol, setup=True)
@@ -192,6 +194,7 @@
         (gmres, 'GMRES'),
         (bicgstab, 'BICGSTAB')]:
     if getattr(d, 'do'+label):
+        assert np.all(np.isfinite(rhs.toarray())), rhs.toarray()
         solver.setPreconditioner(acc)
         solver.setInitialGuess()
         with d.timer('Solve '+label):
@@ -203,14 +206,17 @@
         its.add('Number of iterations '+label, numIter)
         res.add('Residual norm '+label, resNorm)
         resHist.add(label, residuals)
+        print(label, r0, residuals, numIter, resNorm)
     if getattr(d, 'doP'+label):
+        assert np.all(np.isfinite(rhs.toarray())), rhs.toarray()
         solver.setPreconditioner(ml.asPreconditioner(cycle=V), False)
         solver.setInitialGuess()
         with d.timer('Solve P'+label):
             numIter = solver(rhs, x)
         residuals = solver.residuals
         A.residual_py(x, rhs, r)
         resNorm = r.norm(False)
+        print(label, r0, residuals, numIter, resNorm)
         rate.add('Rate of convergence P'+label, (resNorm/r0)**(1/numIter), tested=False if label == 'BICGSTAB' else None)
         its.add('Number of iterations P'+label, numIter, aTol=2 if label == 'BICGSTAB' else None)
         res.add('Residual norm P'+label, resNorm)

multilevelSolver/PyNucleus_multilevelSolver/coarseSolvers_{SCALAR}.pxi

@@ -108,6 +108,7 @@ cdef class {SCALAR_label}coarseSolver({SCALAR_label_lc_}iterative_solver):
     cdef BOOL_t solve_cg(self):
         self.rhs[:] = 0.
         self.overlapsCoarse.receive{SCALAR_label}(self.rhs)
+        print('receive rhs', np.linalg.norm(np.asarray(self.rhs)))
         if self.intraLevelCoarse is not None:
             self.intraLevelCoarse.distribute{SCALAR_label}(self.rhs)
         with self.PLogger.Timer('solveTimeLocal'):
@@ -134,6 +135,7 @@ cdef class {SCALAR_label}coarseSolver({SCALAR_label_lc_}iterative_solver):
         if (self.overlapsFine is not None) or (b.shape[0] != self.Ainv.num_rows):
             with self.PLogger.Timer('solveTime'):
                 if self.inSubdomain and self.canWriteRHS():
+                    print('send rhs', np.linalg.norm(np.asarray(b)))
                     self.sendRHS(b)
                 if self.inCG:
                     ret = self.solve_cg()

multilevelSolver/PyNucleus_multilevelSolver/multigrid_{SCALAR}.pxi

@@ -242,6 +242,8 @@ cdef class {SCALAR_label}multigrid({SCALAR_label_lc_}iterative_solver):
             str label = str(lvlNo)
             FakeTimer Timer = self.PLogger.Timer(label, manualDataEntry=True)
             {SCALAR_label}levelMemory lvl, lvlCoarse
+        print('RHS on ', lvlNo, np.linalg.norm(np.asarray(b)))
+        print('X pre on ', lvlNo, np.linalg.norm(np.asarray(x)))
         if lvlNo == 0:
             Timer.start()
             if isinstance(self.coarse_solver, {SCALAR_label_lc_}iterative_solver):
@@ -289,6 +291,7 @@ cdef class {SCALAR_label}multigrid({SCALAR_label_lc_}iterative_solver):
             lvl.smoother.eval(b, x, postsmoother=True)
             Timer.end()
         Timer.enterData()
+        print('X post on ', lvlNo, np.linalg.norm(np.asarray(x)))
 
     def asPreconditioner(self, INDEX_t maxIter=1, CycleType cycle=V):
         return {SCALAR_label}multigridPreconditioner(self, cycle, maxIter)