Thanks Barry, Next issue:
Consider the following problem:
psiA = cA - 2*cB (1a)
psiC = cC + 2*cB (1b)
cC^2 = k*cA^2*cB (1c)
where psiA, psiC, and k are given. If I solve these five problems using the
standard Newton method:
Case 0:
psiA: 9.400282e-01
psiC: 6.045961e-09
k: 2.000000
cA: 9.400282e-01
cB: 1.555428e-16
cC: 6.045961e-09
Case 1:
psiA: 8.472901e-01
psiC: 7.425271e-09
k: 2.000000
cA: 8.472901e-01
cB: 2.602870e-16
cC: 7.425270e-09
Case 2:
psiA: 8.337199e-01
psiC: 7.831614e-09
k: 2.000000
cA: 8.337199e-01
cB: 2.942675e-16
cC: 7.831613e-09
Case 3:
psiA: 8.268140e-01
psiC: 7.912911e-09
k: 2.000000
cA: 8.268140e-01
cB: 3.029178e-16
cC: 7.912910e-09
Case 4:
psiA: 9.852477e-01
psiC: 7.992223e-09
k: 2.000000
cA: 9.852477e-01
cB: 2.593282e-16
cC: 7.992222e-09
These solutions are "correct", that is if I plug the c's back into
equations (1a)-(1b), i get the original psi's.
Now suppose I use the Variational Inequality such that cA/B/C >= 0, I would
expect to get the exact same solution (since the c's would be non-negative
regardless). But instead I get these solutions:
Case 0:
psiA: 9.400282e-01
psiC: 6.045961e-09
k: 2.000000
cA: 1.318866e-16
cB: 3.097570e-08
cC: 6.045961e-09
Case 1:
psiA: 8.472901e-01
psiC: 7.425271e-09
k: 2.000000
cA: 4.624944e-17
cB: 3.015792e-08
cC: 7.425271e-09
Case 2:
psiA: 8.337199e-01
psiC: 7.831614e-09
k: 2.000000
cA: 1.733276e-16
cB: 3.079856e-08
cC: 7.831614e-09
Case 3:
psiA: 8.268140e-01
psiC: 7.912911e-09
k: 2.000000
cA: 1.721078e-16
cB: 3.061785e-08
cC: 7.912911e-09
Case 4:
psiA: 9.852477e-01
psiC: 7.992223e-09
k: 2.000000
cA: 2.666770e-16
cB: 4.610822e-08
cC: 7.992223e-09
Basically, the cA's shoot down to zero and the cB's are slightly greater
than zero. When I plug the c's into equations (1a) - (1b), I do not get the
correct solution at all. I would expect discrepancies if my c's were
originally negative, but for these five problems, shouldn't VI give the
same answer as the Newton's method?
Attached is the petsc4py code, the two input files containing psiA and
psiC, and the outputs from '-snes_monitor -snes_view
-ksp_monitor_true_residual'. Run as:
python testVI.py testA testC 2 <0/1>
where 0 is Newton, 1 is VI.
Know what's going on here?
Thanks,
Justin
On Wed, Dec 2, 2015 at 1:36 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:
>
> > On Dec 2, 2015, at 1:56 PM, Justin Chang <jychan...@gmail.com> wrote:
> >
> > Barry,
> >
> > So I do not believe there is any problem with the ISEqual() per se,
> because what I am doing is I am solving 5151 different VI problem using the
> same SNES/KSP/PC/Mat. That is, I am not "resetting" these data structures
> once I am done with one problem and move on to the next. If I ran each case
> individually, I get no error; the error comes when I run these problems
> sequentially without resetting the SNES after each problem.
> >
> > I haven't run the C debugger or done any of that yet, but could this be
> a plausible explanation for my error?
>
> This is exactly the issue. Independent of VIs if you change the size of
> a problem you pass to SNES you need to call SNESReset() between each call
> of a different sizes.
>
> > Originally I was thinking about creating/destroying SNES for each
> problem but I was wondering if that might affect performance.
>
> Using SNESReset() is the way you should do it. Creating and destroying
> it each time should only be a tiny, immeasurably slower.
>
> Barry
>
> >
> > Thanks,
> > Justin
> >
> > On Tue, Dec 1, 2015 at 5:58 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:
> >
> > > On Dec 1, 2015, at 5:13 PM, Justin Chang <jychan...@gmail.com> wrote:
> > >
> > > Hi all,
> > >
> > > So I am running into some issues with the SNESVINEWTONRSLS solver. I
> originally had this done in firedrake, but have ported it to petsc4py.
> Attached is the code as well as the two required input files.
> > >
> > > First issue, when I run the program like this:
> > >
> > > python testVI.py psiA_1 psiB_1 2 1
> > >
> > > I get this error:
> > >
> > > Traceback (most recent call last):
> > > File "testVI.py", line 103, in <module>
> > > snes.solve(None,X)
> > > File "PETSc/SNES.pyx", line 520, in petsc4py.PETSc.SNES.solve
> (src/petsc4py.PETSc.c:169677)
> > > petsc4py.PETSc.Error: error code 60
> > > [0] SNESSolve() line 3992 in
> /Users/justin/Software/petsc/src/snes/interface/snes.c
> > > [0] SNESSolve_VINEWTONRSLS() line 507 in
> /Users/justin/Software/petsc/src/snes/impls/vi/rs/virs.c
> > > [0] KSPSetOperators() line 544 in
> /Users/justin/Software/petsc/src/ksp/ksp/interface/itcreate.c
> > > [0] PCSetOperators() line 1170 in
> /Users/justin/Software/petsc/src/ksp/pc/interface/precon.c
> > > [0] Nonconforming object sizes
> > > [0] Cannot change local size of Amat after use old sizes 2 2 new sizes
> 3 3
> > >
> > > No such error occurred when this was ran in firedrake, though I did
> notice that -snes_view did indicate that some of the iterations had 2x2
> instead of 3x3 matrices. Why is this happening? I thought I should be
> solving a 3x3 system each time so why would a 2x2 ever arise?
> >
> > Justin,
> >
> > For VI's the rs (reduced space solver) solves the linearized problem
> on a subset of the variables, thus the size of the linear system may change
> from iteration of Newton to the next iteration of Newton.
> >
> > In the rs solver we have
> >
> > ierr =
> ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal);CHKERRQ(ierr);
> > if (!isequal) {
> > ierr =
> SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr);
> > }
> >
> > so what __should__ happen is that if the size of the reduced space
> changes KSPReset() is called on the KSP allowing the KSP/PC to handle a
> different sized system then before. But why doesn't it work? Why isn't
> KSPReset() called? Somehow the ISEqual is not working.
> >
> > Can you run the C debugger and put a breakpoint in the ISEqual() then
> check the inputs and see if it is correctly setting the isequal to false
> when it should? Each time the vi->IS_inact changes the KSPReset() should
> get called. You can check this in the debugger.
> >
> > Barry
> >
> > >
> > > More issues to come :)
> > >
> > > Thanks,
> > > Justin
> > > <testVI.py><psiA_1><psiB_1>
> >
> >
>
>
#==============================================
# I am solving this chemical reaction system:
#
# 2A - 2C <=> -B
#
# The algebraic formulation:
#
# psiA = cA - 2*cB
# psiC = cC + 2*cB
# cC^2 = keq*cA^2*cB
#
# where psiA and psiC are the totals obtained
# from the advection-diffusion equation. Let
# cA and cC be the primary species and cB be
# the secondary species
#
# HOW TO RUN:
#
# python testVI.py <psiA file> <psiC file> <k> <VI>
#
# <psiA file> - file containing list of psiA
# <psiC file> - file containing list of psiC
# <k> - scalar value for equilibrium constant
# <VI> - Use Variational Inequality (1) or not (0)
#
#==============================================
import sys
#========================
# Initialize
#========================
import petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc
import numpy as np
from math import pow
psiA = np.loadtxt(str(sys.argv[1]))
psiC = np.loadtxt(str(sys.argv[2]))
#psiA = float(sys.argv[1])
#psiC = float(sys.argv[2])
#TOL = 1e-11
k1 = float(sys.argv[3])
opt_VI = int(sys.argv[4])
if psiA.size != psiC.size:
print "Error: mismatch in shape sizes, psiA = " + \
str(psiA.size) + ", psiC = " + str(psiC.size)
exit()
#========================
# Variational inequality
#========================
class Reaction(object):
def __init__(self):
self.psiA = 0
self.psiC = 0
self.k1 = 0
def formData(self, x, a, c, k):
X = x.getArray()
X[0] = 0
X[1] = 0
X[2] = 0
self.psiA = a
self.psiC = c
self.k1 = k
def formFunction(self, snes, x, f):
X = x.getArray(readonly=True)
F = f.getArray()
F[0] = X[0] - 2*X[1] - self.psiA
F[1] = X[2] + 2*X[1] - self.psiC
F[2] = pow(X[2],2) - self.k1*pow(X[0],2)*X[1]
def formJacobian(self, snes, x, J, P):
X = x.getArray(readonly=True)
P.zeroEntries()
P.setValues([0,1,2],[0,1,2],[1,-2,0,0,2,1,-2*self.k1*X[0]*X[1],-self.k1*pow(X[0],2),2*X[2]])
P.assemble()
if J != P:
J.assemble()
return PETSc.Mat.Structure.SAME_NONZERO_PATTERN
pde = Reaction()
snes = PETSc.SNES().create(comm=PETSc.COMM_SELF)
A = PETSc.Mat().createDense(size=[3,3],comm=PETSc.COMM_SELF)
A.setUp()
J = A
F,X = A.createVecs()
snes.setFunction(pde.formFunction,F)
snes.setJacobian(pde.formJacobian,J,J)
if opt_VI == 1:
vilb = X.duplicate()
viub = X.duplicate()
vilb.set(0)
viub.set(10)
snes.setType('vinewtonrsls')
snes.setVariableBounds(vilb,viub)
snes.setFromOptions()
#================================
# Iterate through psi list
#================================
for i in range(0,psiA.size):
pde.formData(X,psiA[i],psiC[i],k1)
if opt_VI == 1:
snes.setVariableBounds(vilb,viub)
snes.solve(None,X)
#X.chop(1e-11)
snes.reset()
x = X.getArray()
print "Case %d:" % i
print "\tpsiA: %e" % psiA[i]
print "\tpsiC: %e" % psiC[i]
print "\tk: %f" % k1
print "\tcA: %e" % x[0]
print "\tcB: %e" % x[1]
print "\tcC: %e" % x[2]
testA
Description: Binary data
testC
Description: Binary data
newton_output
Description: Binary data
vi_output
Description: Binary data
