Since I can't send binary files, I hardcoded the values on main.py. You
only need main.py file to run the code.
Thanks,
On Mon, Nov 12, 2018 at 5:17 PM HeeHo Park <heeho.p...@gmail.com> wrote:
> Hi PETSc,
>
> These Python looks complicated but it really isn't. please run main.py
>
> I am writing up a script that would read petsc binary files extracted from
> one of PFLOTRAN newton iteration step and solve with petsc4py to
> investigate linear solvers.
>
> But I am having trouble solving a matrix that converges PFLOTRAN, let
> alone running a trimmed size (15x15) matrix of A. You'll see in the
> script.
>
> I don't think there is any problem with loading sparse.csr_matrix onto
> PETSc and same with array as I tested with very simple matrix which it
> solves.
>
> Can you figure out why I'm getting inf solutions?
>
> If I use the sparse matrix package to solve it, it solves it very quickly.
>
> x = sla.spsolve(A,b)
>
> Thanks,
>
> --
> HeeHo Daniel Park
>
--
HeeHo Daniel Park
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as sla
import numpy.linalg as nla
#import solve
#import load_data as load
#reload(load)
#reload(solve)
import petsc4py
import sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
#data_path = '/home/heeho/OneDrive/Documents/working/LS_breakdown/simple/no_CheckUpdatePre/'
#data_path = '/home/heeho/OneDrive/Documents/working/LS_breakdown/simple/bjacobi_bicgstab_1core/'
#matrix_name = 'Gjacobian_ts2_tc0_ni0.bin'
#vector_name = 'Gresidual_ts2_tc0_ni0.bin'
#solution_name0 = 'Gxx_ts2_tc0_ni0.bin'
#solution_name1 = 'Gxx_ts2_tc0_ni1.bin'
#A = load.petsc_binary_as_csr(data_path, matrix_name)
#b = load.petsc_binary_vector(data_path, vector_name)
#xx0 = load.petsc_binary_vector(data_path, solution_name0)
#xx1 = load.petsc_binary_vector(data_path, solution_name1)
#A = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]])
#A = sp.csr_matrix(A)
#b = np.random.rand(3)
#A = A[0:15,0:15]
#b = b[0:15]
A = np.array([[ 1.95547767e-14, -3.95323131e-09, 5.40745577e-13,
-3.25907984e-15, 5.54664403e-10, -1.80619170e-15,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 1.95547962e-20, 3.95323131e-09, 5.40746581e-19,
-3.25908310e-21, -5.54664403e-10, -1.80618999e-21,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 4.16896698e-14, 0.00000000e+00, 3.24020244e-05,
-6.94823603e-15, 0.00000000e+00, -5.39999999e-06,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -3.25908610e-15, 5.67439405e-10, -2.81591561e-13,
1.62956619e-14, -3.38578322e-09, 2.57309015e-13,
-3.25907356e-15, 5.54663441e-10, -1.80664029e-15,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -3.25908936e-21, -5.67439405e-10, -2.81591840e-19,
1.62956781e-20, 3.38578324e-09, 2.57309268e-19,
-3.25907681e-21, -5.54663441e-10, -1.80664061e-21,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -6.94846016e-15, 0.00000000e+00, -5.40000159e-06,
3.47409375e-14, 0.00000000e+00, 2.70020226e-05,
-6.94801192e-15, 0.00000000e+00, -5.40000001e-06,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-3.25907986e-15, 5.67441855e-10, -2.81667313e-13,
1.62956330e-14, -3.41133674e-09, 8.24339218e-13,
-3.25906727e-15, 5.54662480e-10, -1.80726833e-15,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-3.25908311e-21, -5.67441855e-10, -2.81667597e-19,
1.62956493e-20, 3.41133674e-09, 8.24339879e-19,
-3.25907053e-21, -5.54662480e-10, -1.80726807e-21,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-6.94823619e-15, 0.00000000e+00, -5.40000154e-06,
3.47406657e-14, 0.00000000e+00, 2.70020258e-05,
-6.94778776e-15, 0.00000000e+00, -5.40000003e-06,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-3.25907357e-15, 5.67445439e-10, -2.81768074e-13,
1.62956042e-14, -3.43690046e-09, 1.39159613e-12,
-3.25906095e-15, 5.54661518e-10, -1.80781860e-15],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-3.25907683e-21, -5.67445439e-10, -2.81768356e-19,
1.62956204e-20, 3.43690046e-09, 1.39159767e-18,
-3.25906421e-21, -5.54661518e-10, -1.80781567e-21],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-6.94801210e-15, 0.00000000e+00, -5.40000156e-06,
3.47403935e-14, 0.00000000e+00, 2.70020290e-05,
-6.94756344e-15, 0.00000000e+00, -5.39999998e-06],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-3.25906729e-15, 5.67448452e-10, -2.81856226e-13,
1.62955751e-14, -3.46247157e-09, 1.95901962e-12],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-3.25907055e-21, -5.67448452e-10, -2.81856509e-19,
1.62955914e-20, 3.46247157e-09, 1.95902204e-18],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-6.94778793e-15, 0.00000000e+00, -5.40000158e-06,
3.47401209e-14, 0.00000000e+00, 2.70020319e-05]])
A = sp.csr_matrix(A)
b = np.array([ -2.78755730e-18, -2.78758937e-24, 4.65216705e-18,
-8.36354149e-18, -8.36364945e-24, 1.40265833e-17,
-1.39417919e-17, -1.39420232e-23, 2.33606221e-17,
-1.95232366e-17, -1.95236906e-23, 3.27508627e-17,
-2.51080558e-17, -2.51089188e-23, 4.20062326e-17])
petsc_A = PETSc.Mat().createAIJ(size=A.shape,
csr=(A.indptr, A.indices, A.data))
petsc_b = PETSc.Vec().createWithArray(b)
petsc_b.assemblyBegin()
petsc_b.assemblyEnd()
ksp = PETSc.KSP().create()
ksp.setOperators(petsc_A)
ksp.setType('preonly')
pc=ksp.getPC()
pc.setType('lu')
#pc.setFactorSolverType('mumps')
#opts = PETSc.Options()
#opts.setValue('-st_pc_factor_mat_solver_package','mumps')
ksp.setFromOptions()
print 'Solving with: ', ksp.getType()
print 'Preconditioning: ', pc.getType()
x = PETSc.Vec().createWithArray(np.zeros(np.size(b)))
x.assemblyBegin()
x.assemblyEnd()
ksp.solve(petsc_b, x)
'''
# ax = load.pflotran_matrix_plot(data_path, matrix_name)
# x = solve.scipy_direct(A, b)
AA,bb,xx = solve.petsc_direct(A, b)
#x = x*9.1778074003631593e-5
'''
'''
data_paths = '/home/heeho/OneDrive/Documents/working/LS_breakdown/simple/direct_1core/'
solution_name0s = 'Gxx_ts2_tc0_ni0.mat'
solution_name1s = 'Gxx_ts2_tc0_ni1.mat'
xx0s = load.matlab_ascii_vector(data_paths, solution_name0s)
xx1s = load.matlab_ascii_vector(data_paths, solution_name1s)
'''
