Re: [petsc-users] Can't figure out why I get inf solutions.
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 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.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00], [ 1.95547962e-20, 3.95323131e-09, 5.40746581e-19, -3.25908310e-21, -5.54664403e-10, -1.80618999e-21, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00], [ 4.16896698e-14, 0.e+00, 3.24020244e-05, -6.94823603e-15, 0.e+00, -5.3999e-06, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+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.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+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.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00], [ -6.94846016e-15, 0.e+00, -5.4159e-06, 3.47409375e-14, 0.e+00, 2.70020226e-05, -6.94801192e-15, 0.e+00, -5.4001e-06, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00], [ 0.e+00, 0.e+00, 0.e+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.e+00, 0.e+00, 0.e+00], [ 0.e+00, 0.e+00, 0.e+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.e+00, 0.e+00, 0.e+00], [ 0.e+00, 0.e+00, 0.e+00, -6.94823619e-15, 0.e+00, -5.4154e-06, 3.47406657e-14, 0.e+00, 2.70020258e-05, -6.94778776e-15, 0.e+00, -5.4003e-06, 0.e+00, 0.e+00, 0.e+00], [ 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, -3.25907357e-15, 5.67445439e-1
Re: [petsc-users] subprocess (block) tolerance.
Barry and Jed, Thank you for your answers. I think I need to learn more about domain decomposition as I am a bit confused. Is it true that we are using BiCGstab here to solve the system of equations, using Additive Schwartz as a domain decomposition preconditioner, and that precondition matrix for each block is preconditioned by ilu[0]? Thanks, On Thu, Oct 4, 2018 at 5:53 PM Smith, Barry F. wrote: > > Since > > KSP Object: (flow_sub_) 1 MPI processes > type: preonly > > this means only a single iteration of the inner solver is used so the > numbers in red are not used. > > You could do something like -flow_ksp_type fgmres -flow_sub_ksp_type gmres > -flow_sub_ksp_rtol 1.e-2 but it wouldn't help maters. Likely the current > values are the best. > >Barry > > > > On Oct 4, 2018, at 5:44 PM, HeeHo Park wrote: > > > > Hi, I'm running PFLOTRAN and in PFLOTRAN, we have flow_ and flow_sub_ > processes. I was wondering what the red underlined values meant (each block > tolerance?) and how to change them (would it affect convergence?). Blue > marked bold values are changed from the default values for linear solvers. > > > > FLOW Linear Solver > >solver: bcgs > >preconditioner: asm > > atol: 1.00E-10 > > rtol: 1.00E-10 > > dtol: 1.00E+04 > > maximum iteration: 1 > > KSP Object: (flow_) 8 MPI processes > > type: bcgs > > maximum iterations=1, initial guess is zero > > tolerances: relative=1e-10, absolute=1e-10, divergence=1. > > left preconditioning > > using PRECONDITIONED norm type for convergence test > > PC Object: (flow_) 8 MPI processes > > type: asm > > Additive Schwarz: total subdomain blocks = 8, amount of overlap = 1 > > Additive Schwarz: restriction/interpolation type - RESTRICT > > [0] number of local blocks = 1 > > [1] number of local blocks = 1 > > [2] number of local blocks = 1 > > [3] number of local blocks = 1 > > [4] number of local blocks = 1 > > [5] number of local blocks = 1 > > [6] number of local blocks = 1 > > [7] number of local blocks = 1 > > Local solve info for each block is in the following KSP and PC > objects: > > - - - - - - - - - - - - - - - - - - > > [0] local block number 0, size = 1389 > > KSP Object: (flow_sub_) 1 MPI processes > > type: preonly > > maximum iterations=1, initial guess is zero > > >>> tolerances: relative=1e-05, absolute=1e-50, divergence=1. > > left preconditioning > > using DEFAULT norm type for convergence test > > PC Object: (flow_sub_) 1 MPI processes > > type: ilu > > PC has not been set up so information may be incomplete > > out-of-place factorization > > 0 levels of fill > > tolerance for zero pivot 2.22045e-14 > > using diagonal shift on blocks to prevent zero pivot [INBLOCKS] > > matrix ordering: natural > > linear system matrix = precond matrix: > > Mat Object: (flow_) 1 MPI processes > > type: seqbaij > > rows=1389, cols=1389, bs=3 > > total: nonzeros=20025, allocated nonzeros=20025 > > total number of mallocs used during MatSetValues calls =0 > > block size is 3 > > - - - - - - - - - - - - - - - - - - > > > > -- > > HeeHo Daniel Park > > -- HeeHo Daniel Park
[petsc-users] subprocess (block) tolerance.
Hi, I'm running PFLOTRAN and in PFLOTRAN, we have flow_ and flow_sub_ processes. I was wondering what the red underlined values meant (each block tolerance?) and how to change them (would it affect convergence?). Blue marked bold values are changed from the default values for linear solvers. FLOW Linear Solver solver: bcgs preconditioner: asm *atol: 1.00E-10* * rtol: 1.00E-10* dtol: 1.00E+04 maximum iteration: 1 KSP Object: (flow_) 8 MPI processes type: bcgs maximum iterations=1, initial guess is zero tolerances: * relative=1e-10, absolute=1e-10*, divergence=1. left preconditioning using PRECONDITIONED norm type for convergence test PC Object: (flow_) 8 MPI processes type: asm Additive Schwarz: total subdomain blocks = 8, amount of overlap = 1 Additive Schwarz: restriction/interpolation type - RESTRICT [0] number of local blocks = 1 [1] number of local blocks = 1 [2] number of local blocks = 1 [3] number of local blocks = 1 [4] number of local blocks = 1 [5] number of local blocks = 1 [6] number of local blocks = 1 [7] number of local blocks = 1 Local solve info for each block is in the following KSP and PC objects: - - - - - - - - - - - - - - - - - - [0] local block number 0, size = 1389 KSP Object: (flow_sub_) 1 MPI processes type: preonly maximum iterations=1, initial guess is zero >>> tolerances: *relative=1e-05, absolute=1e-50*, divergence=1. left preconditioning using DEFAULT norm type for convergence test PC Object: (flow_sub_) 1 MPI processes type: ilu PC has not been set up so information may be incomplete out-of-place factorization 0 levels of fill tolerance for zero pivot 2.22045e-14 using diagonal shift on blocks to prevent zero pivot [INBLOCKS] matrix ordering: natural linear system matrix = precond matrix: Mat Object: (flow_) 1 MPI processes type: seqbaij rows=1389, cols=1389, bs=3 total: nonzeros=20025, allocated nonzeros=20025 total number of mallocs used during MatSetValues calls =0 block size is 3 - - - - - - - - - - - - - - - - - - -- HeeHo Daniel Park
Re: [petsc-users] Fwd: what is the equivalent DMDAVecRestoreArray() function in petsc4py?
Yes, this works. u = da.getVecArray(U) for i in range(mstart, mend): u[i] = np.sin(np.pi*i*6.*h) + 3.*np.sin(np.pi*i*2.*h) The code above also worked without restoreVecArray. I guess the u just points at the array U. I think your code is clearer to understand what is happening. Thank you, On Wed, Feb 14, 2018 at 5:57 PM, Matthew Knepley <knep...@gmail.com> wrote: > On Wed, Feb 14, 2018 at 6:05 PM, HeeHo Park <heeho.p...@gmail.com> wrote: > >> I just found a user group on PETSc website. Can someone please answer the >> question below? >> > > I think it will work using > > with da.getVecArray(U) as u > for i in range(mstart, mend): > u[i] = np.sin(np.pi*i*6.*h) + 3.*np.sin(np.pi*i*2.*h) > > Does it? > > Thanks, > > Matt > > Thanks! >> >> -- Forwarded message -- >> From: HeeHo Park <heeho.p...@gmail.com> >> Date: Wed, Feb 14, 2018 at 5:04 PM >> Subject: what is the equivalent DMDAVecRestoreArray() function in >> petsc4py? >> To: dalc...@gmail.com >> >> >> Hi Lisandro, >> >> I cannot find DMDAVecRestoreArray() equivalent in petsc4py. >> I'm trying to set a 1D initial condition like this. >> >> def initial_conditions(ts, U, appctx): >> da = ts.getDM() >> mstart,xm = da.getCorners() >> mstart = mstart[0] >> xm = xm[0] >> M = da.getSizes()[0] >> h = 1.0/M >> mend = mstart + xm >> >> u = da.getVecArray(U) >> for i in range(mstart, mend): >> u[i] = np.sin(np.pi*i*6.*h) + 3.*np.sin(np.pi*i*2.*h) >> >> da.getVecRestoreArray(u) >> >> Also, is there a better way to ask questions about petsc4py? a forum? or >> google-group? >> >> Thanks, >> >> -- >> HeeHo Daniel Park >> >> >> >> -- >> HeeHo Daniel Park >> > > > > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener > > https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/> > -- HeeHo Daniel Park
[petsc-users] Fwd: what is the equivalent DMDAVecRestoreArray() function in petsc4py?
I just found a user group on PETSc website. Can someone please answer the question below? Thanks! -- Forwarded message -- From: HeeHo Park <heeho.p...@gmail.com> Date: Wed, Feb 14, 2018 at 5:04 PM Subject: what is the equivalent DMDAVecRestoreArray() function in petsc4py? To: dalc...@gmail.com Hi Lisandro, I cannot find DMDAVecRestoreArray() equivalent in petsc4py. I'm trying to set a 1D initial condition like this. def initial_conditions(ts, U, appctx): da = ts.getDM() mstart,xm = da.getCorners() mstart = mstart[0] xm = xm[0] M = da.getSizes()[0] h = 1.0/M mend = mstart + xm u = da.getVecArray(U) for i in range(mstart, mend): u[i] = np.sin(np.pi*i*6.*h) + 3.*np.sin(np.pi*i*2.*h) da.getVecRestoreArray(u) Also, is there a better way to ask questions about petsc4py? a forum? or google-group? Thanks, -- HeeHo Daniel Park -- HeeHo Daniel Park
