Commit dc908316 authored by sfritschi's avatar sfritschi
Browse files

Added debug functionality

parent 5b52dc8e
......@@ -155,9 +155,7 @@ def __amg_pressure_system(pore_dict: Dict[Pore,int], A: Dict, b: Dict, \
bn = data
# iterative solver (tol is relative residual r[k]/r[0] tolerance)
import pyamg
from time import perf_counter
start = perf_counter()
ml = pyamg.ruge_stuben_solver(An) # construct multigrid hierarchy
if resinfo:
rl = list(); j = 0 # residual list, iteration counter
......@@ -168,23 +166,12 @@ def __amg_pressure_system(pore_dict: Dict[Pore,int], A: Dict, b: Dict, \
callback=printres, residuals=rl)
else:
p = ml.solve(bn.T, tol=tol, maxiter=maxiter, cycle='W')
end = perf_counter()
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
# DEBUGGING
if rank == 0:
print("Solution norm pyamg: %.9e" % np.linalg.norm(p))
#print(ml)
print("Time AMG: %e s" % (end - start))
return {pore:p[k] for pore, k in pore_dict.items()}
def __petsc_pressure_system(pore_dict: Dict[Pore,int], A: Dict, b: Dict, \
tol: float, maxiter: int, resinfo: bool = False) -> Dict[Pore,float]:
"""Solve pressure system with PETSc iterative solver in C."""
from time import perf_counter
from mpi4py import MPI
comm = MPI.COMM_WORLD
......@@ -220,13 +207,8 @@ def __petsc_pressure_system(pore_dict: Dict[Pore,int], A: Dict, b: Dict, \
print("Missing dynamic C library; Run 'make all' first")
comm.Abort(1)
start = perf_counter()
p = solve_py(An, bn, tol, maxiter, n)
end = perf_counter()
# DEBUGGING
if rank == root:
print("Solution norm PETSC: %.9e" % np.linalg.norm(p))
print("Time PETSC: %e s" % (end - start))
return {pore:p[k] for pore, k in pore_dict.items()}
def __cg_pressure_system(pore_dict: Dict[Pore,int], A: Dict, b: Dict, \
......@@ -255,14 +237,6 @@ def __cg_pressure_system(pore_dict: Dict[Pore,int], A: Dict, b: Dict, \
p = linalg.cg(An, bn, tol=tol, maxiter=maxiter, callback=printres)[0]
else:
p = linalg.cg(An, bn, tol=tol, maxiter=maxiter)[0]
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
# DEBUGGING
if rank == 0:
print("Solution norm CG: %.9e" % np.linalg.norm(p))
return {pore:p[k] for pore, k in pore_dict.items()}
......@@ -298,14 +272,6 @@ def __ilu_gmres_pressure_system(pore_dict: Dict[Pore,int], A: Dict, b: Dict, \
callback=printres)[0]
else:
p = linalg.gmres(An, bn, M=M, tol=tol, maxiter=maxiter)[0]
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
# DEBUGGING
if rank == 0:
print("Solution norm ILU: %.9e" % np.linalg.norm(p))
return {pore:p[k] for pore, k in pore_dict.items()}
......
......@@ -9,6 +9,9 @@ export PETSC_LIB
all: setup.py
python3 $< build_ext -i
debug: setup.py
python3 $< build_ext -i --debug
run:
mpirun -n 4 python3 main.py
......
......@@ -34,6 +34,12 @@ link_args += os.popen(mpi + " --showme:link").read().strip().split(" ")
link_args += PETSC_LIB.strip().split(" ")
sources = ["cwrapper/wrapper.pyx", "src/solve.c"]
# Set debug macro
macros = [('DEBUG', '0')]
if "--debug" in sys.argv:
macros = [('DEBUG', '1')]
# Remove argument again
sys.argv.remove("--debug")
setup(
name = "wrapper",
......@@ -41,7 +47,8 @@ setup(
[Extension("wrapper",
sources=sources,
extra_compile_args=compile_args,
extra_link_args=link_args)], \
extra_link_args=link_args,
define_macros=macros)], \
compiler_directives={'language_level': 3},
)
)
......@@ -3,6 +3,14 @@
#include <stdio.h> // printf
#include <stdlib.h> // malloc
// Ignore macros in release build
#ifndef DEBUG
#undef CHKERRQ
#define CHKERRQ(ierr)
#undef CHKERRMPI
#define CHKERRMPI(ierr)
#endif
// Call this before running the function
PetscErrorCode init() {
PetscErrorCode ierr;
......@@ -94,14 +102,14 @@ PetscErrorCode solve(const PetscScalar *data, const PetscInt *col_indices,
ierr = KSPSolve(ksp, bn, x); CHKERRQ(ierr);
// DEBUGGING: Monitor iteration number/convergence
PetscInt its;
ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "Iterations %d\n", its);
const char *convergence_reason;
ierr = KSPGetConvergedReasonString(ksp, &convergence_reason); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "Reason: %s\n", convergence_reason);
if (DEBUG) {
PetscInt its;
ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "Iterations %d\n", its);
const char *convergence_reason;
ierr = KSPGetConvergedReasonString(ksp, &convergence_reason); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "Reason: %s\n", convergence_reason);
}
// Displacements in global array (inscan of n_locals array)
PetscMPIInt *displs = NULL;
......
......@@ -23,7 +23,7 @@ def main():
inpores.add(pore)
elif pore.label[:3] == 'out':
outpores.add(pore)
print("Pores: %d In-pores: %d Out-pores: %d" % (len(basenet.pores), \
len(inpores), \
len(outpores)))
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment