diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 38fbabaa..77bb7152 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -107,7 +107,7 @@ jobs: - name: Run tests if: always() - run: python3 -m pytest --junit-xml=test-results-${{ env.BUILD_IDENTIFIER }}.xml tests/ + run: python3 -m pytest --junit-xml=test-results-${{ env.BUILD_IDENTIFIER }}.xml -k "testParallelGMG[4-interval-P3-False]" tests/ - name: Run flake8 if: always() diff --git a/Makefile b/Makefile index d2d512cd..e10795c9 100644 --- a/Makefile +++ b/Makefile @@ -130,7 +130,8 @@ list-tests: .PHONY: tests tests: - $(PYTHON) -m pytest -rA --html=$(TEST_RESULTS) --self-contained-html tests/ + # $(PYTHON) -m pytest -rA --html=$(TEST_RESULTS) --self-contained-html tests/ + $(PYTHON) -m pytest -rA --html=$(TEST_RESULTS) -k "testParallelGMG[4-interval-P3-False]" --self-contained-html tests/ prereq: $(PYTHON) -m pip install $(PIP_FLAGS) $(PIP_INSTALL_FLAGS) setuptools wheel Cython cython numpy scipy matplotlib pyyaml h5py pybind11 MeshPy tabulate modepy mpi4py pyamg meshio diff --git a/base/PyNucleus_base/solvers.pyx b/base/PyNucleus_base/solvers.pyx index f9503771..63c55100 100644 --- a/base/PyNucleus_base/solvers.pyx +++ b/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): diff --git a/drivers/runParallelGMG.py b/drivers/runParallelGMG.py index 3df1b2aa..0539ffc9 100755 --- a/drivers/runParallelGMG.py +++ b/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,7 +206,9 @@ 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): @@ -211,6 +216,7 @@ 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) diff --git a/fem/PyNucleus_fem/algebraicOverlaps.pyx b/fem/PyNucleus_fem/algebraicOverlaps.pyx index 45a9a076..a70fff60 100644 --- a/fem/PyNucleus_fem/algebraicOverlaps.pyx +++ b/fem/PyNucleus_fem/algebraicOverlaps.pyx @@ -1170,6 +1170,7 @@ cdef class algebraicOverlapManager: algebraicOverlapOneSidedGet ov1S algebraicOverlapOneSidedPut ov1SP algebraicOverlapOneSidedPutLockAll ov1SPLA + print('pre-receive', self.comm.rank, np.linalg.norm(np.asarray(return_vec))) requestsOneSidedGet = [] requestsOneSidedPut = [] requestsOneSidedPutLockAll = [] @@ -1215,6 +1216,7 @@ cdef class algebraicOverlapManager: ov1SPLA = self.overlaps[subdomainNo] ov1SPLA.receive(return_vec, vecNo=vecNo, flushOp=flushOp) ov1SPLA.accumulateProcess(return_vec, vecNo=vecNo) + print('post-receive', self.comm.rank, np.linalg.norm(np.asarray(return_vec))) if setBarrier: self.comm.Barrier() diff --git a/multilevelSolver/PyNucleus_multilevelSolver/coarseSolvers_{SCALAR}.pxi b/multilevelSolver/PyNucleus_multilevelSolver/coarseSolvers_{SCALAR}.pxi index e690c6cd..f3ceac6d 100644 --- a/multilevelSolver/PyNucleus_multilevelSolver/coarseSolvers_{SCALAR}.pxi +++ b/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,11 +135,13 @@ 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() if self.inSubdomain: self.getSolution(x) + print('send solution', np.linalg.norm(np.asarray(x))) else: with self.PLogger.Timer('solveTime'): if isinstance(self.Ainv, {SCALAR_label_lc_}iterative_solver): diff --git a/multilevelSolver/PyNucleus_multilevelSolver/multigrid_{SCALAR}.pxi b/multilevelSolver/PyNucleus_multilevelSolver/multigrid_{SCALAR}.pxi index 9fa76948..b8ccf628 100644 --- a/multilevelSolver/PyNucleus_multilevelSolver/multigrid_{SCALAR}.pxi +++ b/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, 'rank', self.overlap.comm.rank, np.linalg.norm(np.asarray(b))) + print('X pre on ', lvlNo, 'rank', self.overlap.comm.rank, 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, 'rank', self.overlap.comm.rank, np.linalg.norm(np.asarray(x))) def asPreconditioner(self, INDEX_t maxIter=1, CycleType cycle=V): return {SCALAR_label}multigridPreconditioner(self, cycle, maxIter)