Ozzy, Thanks for the update. I have changed the checking of the f() as suggested by the cited reference in master.
Barry > On Jun 1, 2015, at 4:07 AM, Asbjørn Nilsen Riseth <[email protected]> > wrote: > > Hi again Barry, > > I sorted out the jacobian issue, and made a comparison between the two > definitions of the active set. > The active set with strict inequality takes the same or fewer Newton steps > than the current petsc code. With a larger search space, I expect this to > happen. snes_vi_monitor logs comparing the two are shown below. > > I could submit a pull request with the change, but we should probably > consider: > 1) Whether this active set definition is consistent with newtonssls > 2) The linear systems to solve becomes larger, so for some cases the overall > performance might not improve so much. > 3) For more flexibility, we could add an option to decide whether to use a > strict inequality or not. This would sort out 1) and 2). > > I don't have much experience with the petsc codebase though, so adding > options might take me some time. > > > Ozzy > > > __ log using my patch __ > 0 SNES VI Function norm 7.796491981333e+02 Active lower constraints 0/18 > upper constraints 0/0 Percent of total 0 Percent of bounded 0 > 1 SNES VI Function norm 2.405400030748e+02 Active lower constraints 0/16 > upper constraints 0/0 Percent of total 0 Percent of bounded 0 > 2 SNES VI Function norm 2.145739795389e+02 Active lower constraints 0/17 > upper constraints 0/0 Percent of total 0 Percent of bounded 0 > 3 SNES VI Function norm 1.942498283668e+02 Active lower constraints 0/13 > upper constraints 0/0 Percent of total 0 Percent of bounded 0 > 4 SNES VI Function norm 1.834306037299e+01 Active lower constraints 0/11 > upper constraints 0/0 Percent of total 0 Percent of bounded 0 > 5 SNES VI Function norm 1.724597091463e+01 Active lower constraints 0/11 > upper constraints 0/0 Percent of total 0 Percent of bounded 0 > 6 SNES VI Function norm 4.210027533399e-02 Active lower constraints 0/10 > upper constraints 0/0 Percent of total 0 Percent of bounded 0 > 7 SNES VI Function norm 3.014124871281e-07 Active lower constraints 0/10 > upper constraints 0/0 Percent of total 0 Percent of bounded 0 > SNES Object:(firedrake_snes_0_) 1 MPI processes > type: vinewtonrsls > maximum iterations=20, maximum function evaluations=10000 > tolerances: relative=0, absolute=1e-06, solution=0 > total number of linear solver iterations=7 > total number of function evaluations=22 > norm schedule ALWAYS > SNESLineSearch Object: (firedrake_snes_0_) 1 MPI processes > type: l2 > maxstep=1.000000e+08, minlambda=1.000000e-12 > tolerances: relative=1.000000e-08, absolute=1.000000e-15, > lambda=1.000000e-08 > maximum iterations=1 > KSP Object: (firedrake_snes_0_) 1 MPI processes > type: preonly > maximum iterations=10000, initial guess is zero > tolerances: relative=1e-05, absolute=1e-50, divergence=10000 > left preconditioning > using NONE norm type for convergence test > PC Object: (firedrake_snes_0_) 1 MPI processes > type: lu > LU: out-of-place factorization > tolerance for zero pivot 2.22045e-14 > matrix ordering: nd > factor fill ratio given 5, needed 1.54545 > Factored matrix follows: > Mat Object: 1 MPI processes > type: seqaij > rows=36, cols=36 > package used to perform factorization: petsc > total: nonzeros=816, allocated nonzeros=816 > total number of mallocs used during MatSetValues calls =0 > using I-node routines: found 15 nodes, limit used is 5 > linear system matrix = precond matrix: > Mat Object: 1 MPI processes > type: seqaij > rows=36, cols=36 > total: nonzeros=528, allocated nonzeros=528 > total number of mallocs used during MatSetValues calls =0 > not using I-node routines > -------------------------------------------------- > > __ log using the original petsc code __ > 0 SNES VI Function norm 7.796491981333e+02 Active lower constraints 12/18 > upper constraints 0/0 Percent of total 0.333333 Percent of bounded 0.333333 > 1 SNES VI Function norm 2.630718602212e+02 Active lower constraints 12/16 > upper constraints 0/0 Percent of total 0.333333 Percent of bounded 0.333333 > 2 SNES VI Function norm 2.363417090057e+02 Active lower constraints 12/17 > upper constraints 0/0 Percent of total 0.333333 Percent of bounded 0.333333 > 3 SNES VI Function norm 1.902271040685e+01 Active lower constraints 12/14 > upper constraints 0/0 Percent of total 0.333333 Percent of bounded 0.333333 > 4 SNES VI Function norm 1.866193366410e+01 Active lower constraints 12/12 > upper constraints 0/0 Percent of total 0.333333 Percent of bounded 0.333333 > 5 SNES VI Function norm 1.865568900723e+01 Active lower constraints 12/12 > upper constraints 0/0 Percent of total 0.333333 Percent of bounded 0.333333 > 6 SNES VI Function norm 2.182461654877e+02 Active lower constraints 10/12 > upper constraints 0/0 Percent of total 0.277778 Percent of bounded 0.277778 > 7 SNES VI Function norm 2.575010971279e-01 Active lower constraints 10/11 > upper constraints 0/0 Percent of total 0.277778 Percent of bounded 0.277778 > 8 SNES VI Function norm 1.056372578821e-05 Active lower constraints 10/10 > upper constraints 0/0 Percent of total 0.277778 Percent of bounded 0.277778 > 9 SNES VI Function norm 4.368019257866e-11 Active lower constraints 10/10 > upper constraints 0/0 Percent of total 0.277778 Percent of bounded 0.277778 > SNES Object:(firedrake_snes_0_) 1 MPI processes > type: vinewtonrsls > maximum iterations=20, maximum function evaluations=10000 > tolerances: relative=0, absolute=1e-06, solution=0 > total number of linear solver iterations=9 > total number of function evaluations=28 > norm schedule ALWAYS > SNESLineSearch Object: (firedrake_snes_0_) 1 MPI processes > type: l2 > maxstep=1.000000e+08, minlambda=1.000000e-12 > tolerances: relative=1.000000e-08, absolute=1.000000e-15, > lambda=1.000000e-08 > maximum iterations=1 > KSP Object: (firedrake_snes_0_) 1 MPI processes > type: preonly > maximum iterations=10000, initial guess is zero > tolerances: relative=1e-05, absolute=1e-50, divergence=10000 > left preconditioning > using NONE norm type for convergence test > PC Object: (firedrake_snes_0_) 1 MPI processes > type: lu > LU: out-of-place factorization > tolerance for zero pivot 2.22045e-14 > matrix ordering: nd > factor fill ratio given 5, needed 1.57895 > Factored matrix follows: > Mat Object: 1 MPI processes > type: seqaij > rows=26, cols=26 > package used to perform factorization: petsc > total: nonzeros=420, allocated nonzeros=420 > total number of mallocs used during MatSetValues calls =0 > using I-node routines: found 11 nodes, limit used is 5 > linear system matrix = precond matrix: > Mat Object: 1 MPI processes > type: seqaij > rows=26, cols=26 > total: nonzeros=266, allocated nonzeros=266 > total number of mallocs used during MatSetValues calls =0 > not using I-node routines > > > On Sat, 30 May 2015 at 01:07 Asbjørn Nilsen Riseth <[email protected]> > wrote: > Hey Barry, > > thanks for the offer to have a look at the code. I ran SNESTEST today: > user-defined failed, 1.0 and -1.0 seemed to work fine. My first step will > have to be to find out what's wrong with my jacobian. If I've still got > issues after that, I'll try to set up an easy-to-experiment code > > The code is a DG0 FVM formulation set up in Firedrake (a "fork" of FEniCS). I > was assuming UFL would sort the Jacobian for me. > Lesson learnt: always do a SNESTEST. > > > Ozzy > > On Fri, 29 May 2015 at 19:21 Barry Smith <[email protected]> wrote: > > > On May 29, 2015, at 4:19 AM, Asbjørn Nilsen Riseth <[email protected]> > > wrote: > > > > Barry, > > > > I changed the code, but it only makes a difference at the order of 1e-10 - > > so that's not the cause. I've attached (if that's appropriate?) the patch > > in case anyone is interested. > > > > Investigating the function values, I see that the first Newton step goes > > towards the expected solution for my function values. Then it shoots back > > close to the initial conditions. > > When does it "shoot back close to the initial conditions"? At the second > Newton step? If so is the residual norm still smaller at the second step? > > > At the time the solver hits tolerance for the inactive set; the function > > value is 100-1000 at some of the active set indices. > > Reducing the timestep by an order of magnitude shows the same behavior for > > the first two timesteps. > > > > Maybe VI is not the right approach. The company I work with seem to just be > > projecting negative values. > > The VI solver is essentially just a "more sophisticated" version of > "projecting negative values" so should not work worse then an ad hoc method > and should generally work better (sometimes much better). > > Is your code something simple you could email me to play with or is it a > big application that would be hard for me to experiment with? > > > > Barry > > > > > I'll investigate further. > > > > Ozzy > > > > > > On Thu, 28 May 2015 at 20:26 Barry Smith <[email protected]> wrote: > > > > Ozzy, > > > > I cannot say why it is implemented as >= (could be in error). Just > > try changing the PETSc code (just run make gnumake in the PETSc directory > > after you change the source to update the library) and see how it affects > > your code run. > > > > Barry > > > > > On May 28, 2015, at 3:15 AM, Asbjørn Nilsen Riseth > > > <[email protected]> wrote: > > > > > > Dear PETSc developers, > > > > > > Is the active set in NewtonRSLS defined differently from the reference* > > > you give in the documentation on purpose? > > > The reference defines the active set as: > > > x_i = 0 and F_i > 0, > > > whilst the PETSc code defines it as x_i = 0 and F_i >= 0 (vi.c: 356) : > > > !((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || > > > (PetscRealPart(f[i]) < 0.0) > > > So PETSc freezes the variables if f[i] == 0. > > > > > > I've been using the Newton RSLS method to ensure positivity in a > > > subsurface flow problem I'm working on. My solution stays almost constant > > > for two timesteps (seemingly independent of the size of the timestep), > > > before it goes towards the expected solution. > > > From my initial conditions, certain variables are frozen because x_i = 0 > > > and f[i] = 0, and I was wondering if that could be the cause of my issue. > > > > > > > > > *: > > > - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for > > > Large-Scale Applications, Optimization Methods and Software, 21 (2006). > > > > > > > > > Regards, > > > Ozzy > > > > <0001-Define-active-and-inactive-sets-correctly.patch> >
