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. 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. 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 > >
From 951f09745b3b633241d9751138e8d78d5384e59f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Asbj=C3=B8rn=20Nilsen=20Riseth?= <[email protected]> Date: Wed, 27 May 2015 18:51:13 +0100 Subject: [PATCH] Define active and inactive sets correctly --- src/snes/impls/vi/vi.c | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/snes/impls/vi/vi.c b/src/snes/impls/vi/vi.c index c2d4b27..f20c550 100644 --- a/src/snes/impls/vi/vi.c +++ b/src/snes/impls/vi/vi.c @@ -67,14 +67,14 @@ PetscErrorCode SNESVIComputeInactiveSetIS(Vec upper,Vec lower,Vec X,Vec F,IS *in ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); /* Compute inactive set size */ for (i=0; i < nlocal; i++) { - if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) nloc_isact++; + if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) >= 0.0))) nloc_isact++; } ierr = PetscMalloc1(nloc_isact,&idx_act);CHKERRQ(ierr); /* Set inactive set indices */ for (i=0; i < nlocal; i++) { - if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) idx_act[i1++] = ilow+i; + if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) >= 0.0))) idx_act[i1++] = ilow+i; } /* Create inactive set IS */ @@ -137,9 +137,9 @@ PetscErrorCode SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dumm rnorm = 0.0; for (i=0; i<n; i++) { - if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]); - else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + 1.e-8 && PetscRealPart(f[i]) >= 0.0) act[0]++; - else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - 1.e-8 && PetscRealPart(f[i]) <= 0.0) act[1]++; + if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) >= 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]); + else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + 1.e-8 && PetscRealPart(f[i]) > 0.0) act[0]++; + else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - 1.e-8 && PetscRealPart(f[i]) < 0.0) act[1]++; else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Can never get here"); } @@ -346,14 +346,14 @@ PetscErrorCode SNESVIGetActiveSetIS(SNES snes,Vec X,Vec F,IS *ISact) ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); /* Compute active set size */ for (i=0; i < nlocal;i++) { - if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) nloc_isact++; + if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) >= 0.0))) nloc_isact++; } ierr = PetscMalloc1(nloc_isact,&idx_act);CHKERRQ(ierr); /* Set active set indices */ for (i=0; i < nlocal; i++) { - if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) idx_act[i1++] = ilow+i; + if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) >= 0.0))) idx_act[i1++] = ilow+i; } /* Create active set IS */ @@ -397,7 +397,7 @@ PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X, PetscReal *f ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); rnorm = 0.0; for (i=0; i<n; i++) { - if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]); + if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) >= 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]); } ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); ierr = VecRestoreArrayRead(snes->xl,&xl);CHKERRQ(ierr); -- 1.9.1
