Re: [petsc-dev] PetscSF in Fortran

[Adrian Croucher]

hi


On 16/10/17 15:16, Jed Brown wrote:

Adrian Croucher  writes:


So do you think the SF stuff in petsc/finclude/petscis.h should be taken
out, and put into a new petsc/finclude/petscsf.h ?

I think that's desirable for symmetry with include/petscsf.h, but it
isn't important for your contribution.


I have got Fortran bindings for PetscSFGetGraph() and PetscSFSetGraph() 
working. First I tried separating the existing SF stuff out of the IS 
modules, but ran into dependency problems in other modules. I wasn't too 
confident about the best way to resolve those without breaking things, 
so I tried it again just doing minimal changes to the existing SF stuff- 
essentially just adding Fortran wrappers for PetscSFGetGraph() and 
PetscSFSetGraph() into src/vec/f90-mod/petscis.h90- and that seemed to 
work fine. Maybe someone more familiar with the setup could separate out 
the SF stuff from IS sometime.



So, now I'm trying to add Fortran bindings for PetscSFBcastBegin() and 
PetscSFBcastEnd().


From the C side I have added the following into 
src/vec/is/sf/interface/f90-custom/zsff90.c:


PETSC_EXTERN void PETSC_STDCALL petscsfbcastbegin_(PetscSF *sf, 
MPI_Datatype *unit, F90Array1d *rptr, F90Array1d *lptr , int *ierr 
PETSC_F90_2PTR_PROTO(rptrd) PETSC_F90_2PTR_PROTO(lptrd))

{
  PetscDataType ptype;
  const void* rootdata;
  void* leafdata;

  *ierr = PetscMPIDataTypeToPetscDataType(*unit, );if (*ierr) return;
  *ierr = F90Array1dAccess(rptr, ptype, (void**)  
PETSC_F90_2PTR_PARAM(rptrd));if (*ierr) return;
  *ierr = F90Array1dAccess(lptr, ptype, (void**)  
PETSC_F90_2PTR_PARAM(lptrd));if (*ierr) return;


  *ierr = PetscSFBcastBegin(*sf, *unit, rootdata, leafdata);

}

and similarly for petscsfbcastend_(). Does this look plausible?

Then some wrappers need to be added to src/vec/f90-mod/petscis.h90. I am 
not sure how to do those.


The difficulty is in declaring the arrays that are passed in, which can 
be of various types. In C they are declared as void*, but I'm not sure 
what to do with that in Fortran. I can't seem to find any other example 
wrappers in PETSc to model it on either. Any suggestions?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-dev] TSComputeIFunction(): User or PETSc responsible for setting residual vector to zero?

[Brad Aagaard]

The PETSc manual pages do not mention whether or not the user is 
responsible for setting the entries in the function vector to zero in 
the user-defined ifunction and rhsfunction.


In looking at the code it is clear TSComputeIFunction() does not set the 
function vector entries to zero. For the function to be computed 
correctly with the current code, the user must zero the function Vec 
passed in to both the ifunction and rhsfunction. This seems like a bug 
to me, because TSComputeIFunction() could potentially use a single 
global Vec in computing the function from the ifunction and rhsfunction, 
so a user zeroing the entries would result in incorrect values.


It seems to be the code should behave like:

TSComputeIFunction calls VecSet(residual, 0.0) BEFOERE calling ifunction 
and rhsfunction.


User ifunction calls scatter with ADD_VALUES
User rhsfunction calls scatter with ADD_VALUES

This way the user code is agnostic as to whether TSComputeIFunction() 
passes in separate vecs to ifunction and rhsfunction or the same vec 
with the proper -1 scaling done in between.


Whether or not TSComputeIFunction() is changed to zero out the function 
Vec, I think the manual page should be updated to indicate to the user 
whether the user function or PETSc function is responsible for zeroing 
the function Vec.


Thanks,
Brad

Re: [petsc-dev] [petsc-checkbuilds] PETSc blame digest (next) 2017-10-17

[Matthew Knepley]

On Tue, Oct 17, 2017 at 4:57 PM, Satish Balay  wrote:

> pushed to knepley/feature-adaptor-plex
>
> https://bitbucket.org/petsc/petsc/commits/54c8e91ecbda26974e0050fcf2abbb
> a3ecc87a6b
>
> For eg: http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/make_next_arch-cuda-double_es.log
>
> The above fixes the following error that has been aborting nightlybuilds
> for past couple of days.
>

Crap! I did not push the finclude thing. Thanks for fixing it.

I did not know about the bfort thing.

I will move PetscGlobalMinMax() to sys. I wanted to be careful and try to
replace everyplace we do this
in PETSc with this function.

  Thanks,

 Matt


> >>>
> /sandbox/petsc/petsc.clone/include/../src/snes/f90-mod/
> ftn-auto-interfaces/petscsnes.h90:574.7:
> Included at /sandbox/petsc/petsc.clone/src/snes/f90-mod/petscsnesmod.
> F:45:
>
>DMAdaptationStrategy c ! DMAdaptationStrategy
>1
> Error: Unclassifiable statement at (1)
> make[2]: *** [arch-cuda-double/obj/src/snes/f90-mod/petscsnesmod.o] Error
> 1
> <<
>
> Also fixes this warning
>
> 
> /sandbox/petsc/petsc.clone/src/snes/utils/ftn-auto/dmadaptf.c: In
> function ‘petscglobalminmax_’:
> /sandbox/petsc/petsc.clone/src/snes/utils/ftn-auto/dmadaptf.c:125:1:
> warning: implicit declaration of function ‘PetscGlobalMinMax’
> [-Wimplicit-function-declaration]
>  *__ierr = PetscGlobalMinMax(
>  ^
> <<
>
>
> Satish
>
>
> On Tue, 17 Oct 2017, Matthew Knepley wrote:
>
> > I already pushed a fix this morning. What did you push?
> >
> >   Matt
> >
> > On Tue, Oct 17, 2017 at 3:59 PM, Satish Balay  wrote:
> >
> > > This has been broken for a couple of days. I pushed a fix [hoping
> > > builds don't abort anymore] - but there are more things to fix..
> > >
> > > Satish
> > >
> > > On Tue, 17 Oct 2017, PETSc checkBuilds wrote:
> > >
> > > >
> > > >
> > > > Dear PETSc developer,
> > > >
> > > > This email contains listings of contributions attributed to you by
> > > > `git blame` that caused compiler errors or warnings in PETSc
> automated
> > > > testing.  Follow the links to see the full log files. Please attempt
> to
> > > fix
> > > > the issues promptly or let us know at petsc-dev@mcs.anl.gov if you
> are
> > > unable
> > > > to resolve the issues.
> > > >
> > > > Thanks,
> > > >   The PETSc development team
> > > >
> > > > 
> > > >
> > > > warnings attributed to commit https://bitbucket.org/petsc/
> > > petsc/commits/5675c17
> > > > DM: Added DMAdaptor
> > > >
> > > >   src/snes/utils/dmadapt.c:610
> > > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > > 2017/10/17/build_next_arch-linux-c89_el6.log]
> > > >   /sandbox/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:610:5:
> > > error: expected expression before '/' token
> > > >
> > > >   src/snes/utils/dmadapt.c:611
> > > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > > 2017/10/17/build_next_arch-linux-c89_el6.log]
> > > >   /sandbox/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:611:5:
> > > error: expected expression before '/' token
> > > >
> > > >   src/snes/utils/dmadapt.c:710
> > > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > > 2017/10/17/build_next_arch-osx-10.6-cxx-pkgs-opt_ipro.log]
> > > >   /Users/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:710:16:
> > > error: no matching function for call to 'PetscFECreateDefault'
> > > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > > 2017/10/17/build_next_arch-freebsd-cxx-pkgs-opt_wii.log]
> > > >   /usr/home/balay/petsc.clone-2/src/snes/utils/dmadapt.c:710:
> 110:
> > > error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> > > 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> > > const char*, PetscInt, _p_PetscFE**)'
> > > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > > 2017/10/17/build_next_arch-linux-xsdk-dbg_cg.log]
> > > >   /sandbox/petsc/petsc.clone/src/snes/utils/dmadapt.c:710:110:
> > > error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> > > 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> > > const char*, PetscInt, _p_PetscFE**)'
> > > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > > 2017/10/17/build_next_arch-osx-10.6-cxx-cmplx-pkgs-dbg_ipro.log]
> > > >   /Users/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:16:
> > > error: no matching function for call to 'PetscFECreateDefault'
> > > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > > 2017/10/17/build_next_arch-linux-opt-cxx-quad_grind.log]
> > > >   /sandbox/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:110:
> > > error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> > > 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> > > const char*, PetscInt, _p_PetscFE**)'
> > > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > > 

Re: [petsc-dev] [petsc-checkbuilds] PETSc blame digest (next) 2017-10-17

[Satish Balay]

pushed to knepley/feature-adaptor-plex

https://bitbucket.org/petsc/petsc/commits/54c8e91ecbda26974e0050fcf2abbba3ecc87a6b

For eg: 
http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/make_next_arch-cuda-double_es.log

The above fixes the following error that has been aborting nightlybuilds for 
past couple of days.

>>>
/sandbox/petsc/petsc.clone/include/../src/snes/f90-mod/ftn-auto-interfaces/petscsnes.h90:574.7:
Included at /sandbox/petsc/petsc.clone/src/snes/f90-mod/petscsnesmod.F:45:

   DMAdaptationStrategy c ! DMAdaptationStrategy
   1
Error: Unclassifiable statement at (1)
make[2]: *** [arch-cuda-double/obj/src/snes/f90-mod/petscsnesmod.o] Error 1
<<

Also fixes this warning


/sandbox/petsc/petsc.clone/src/snes/utils/ftn-auto/dmadaptf.c: In function 
‘petscglobalminmax_’:
/sandbox/petsc/petsc.clone/src/snes/utils/ftn-auto/dmadaptf.c:125:1: warning: 
implicit declaration of function ‘PetscGlobalMinMax’ 
[-Wimplicit-function-declaration]
 *__ierr = PetscGlobalMinMax(
 ^
<<


Satish


On Tue, 17 Oct 2017, Matthew Knepley wrote:

> I already pushed a fix this morning. What did you push?
> 
>   Matt
> 
> On Tue, Oct 17, 2017 at 3:59 PM, Satish Balay  wrote:
> 
> > This has been broken for a couple of days. I pushed a fix [hoping
> > builds don't abort anymore] - but there are more things to fix..
> >
> > Satish
> >
> > On Tue, 17 Oct 2017, PETSc checkBuilds wrote:
> >
> > >
> > >
> > > Dear PETSc developer,
> > >
> > > This email contains listings of contributions attributed to you by
> > > `git blame` that caused compiler errors or warnings in PETSc automated
> > > testing.  Follow the links to see the full log files. Please attempt to
> > fix
> > > the issues promptly or let us know at petsc-dev@mcs.anl.gov if you are
> > unable
> > > to resolve the issues.
> > >
> > > Thanks,
> > >   The PETSc development team
> > >
> > > 
> > >
> > > warnings attributed to commit https://bitbucket.org/petsc/
> > petsc/commits/5675c17
> > > DM: Added DMAdaptor
> > >
> > >   src/snes/utils/dmadapt.c:610
> > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > 2017/10/17/build_next_arch-linux-c89_el6.log]
> > >   /sandbox/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:610:5:
> > error: expected expression before '/' token
> > >
> > >   src/snes/utils/dmadapt.c:611
> > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > 2017/10/17/build_next_arch-linux-c89_el6.log]
> > >   /sandbox/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:611:5:
> > error: expected expression before '/' token
> > >
> > >   src/snes/utils/dmadapt.c:710
> > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > 2017/10/17/build_next_arch-osx-10.6-cxx-pkgs-opt_ipro.log]
> > >   /Users/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:710:16:
> > error: no matching function for call to 'PetscFECreateDefault'
> > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > 2017/10/17/build_next_arch-freebsd-cxx-pkgs-opt_wii.log]
> > >   /usr/home/balay/petsc.clone-2/src/snes/utils/dmadapt.c:710:110:
> > error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> > 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> > const char*, PetscInt, _p_PetscFE**)'
> > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > 2017/10/17/build_next_arch-linux-xsdk-dbg_cg.log]
> > >   /sandbox/petsc/petsc.clone/src/snes/utils/dmadapt.c:710:110:
> > error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> > 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> > const char*, PetscInt, _p_PetscFE**)'
> > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > 2017/10/17/build_next_arch-osx-10.6-cxx-cmplx-pkgs-dbg_ipro.log]
> > >   /Users/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:16:
> > error: no matching function for call to 'PetscFECreateDefault'
> > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > 2017/10/17/build_next_arch-linux-opt-cxx-quad_grind.log]
> > >   /sandbox/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:110:
> > error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> > 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> > const char*, PetscInt, _p_PetscFE**)'
> > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > 2017/10/17/build_next_arch-freebsd-cxx-cmplx-pkgs-dbg_wii.log]
> > >   /usr/home/balay/petsc.clone/src/snes/utils/dmadapt.c:710:110:
> > error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> > 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> > const char*, PetscInt, _p_PetscFE**)'
> > > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> > 2017/10/17/build_next_arch-linux-pkgs-cxx-mlib_el6.log]
> > >   /sandbox/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:110:
> > error: cannot convert 'bool' to 

Re: [petsc-dev] [petsc-checkbuilds] PETSc blame digest (next) 2017-10-17

[Matthew Knepley]

I already pushed a fix this morning. What did you push?

  Matt

On Tue, Oct 17, 2017 at 3:59 PM, Satish Balay  wrote:

> This has been broken for a couple of days. I pushed a fix [hoping
> builds don't abort anymore] - but there are more things to fix..
>
> Satish
>
> On Tue, 17 Oct 2017, PETSc checkBuilds wrote:
>
> >
> >
> > Dear PETSc developer,
> >
> > This email contains listings of contributions attributed to you by
> > `git blame` that caused compiler errors or warnings in PETSc automated
> > testing.  Follow the links to see the full log files. Please attempt to
> fix
> > the issues promptly or let us know at petsc-dev@mcs.anl.gov if you are
> unable
> > to resolve the issues.
> >
> > Thanks,
> >   The PETSc development team
> >
> > 
> >
> > warnings attributed to commit https://bitbucket.org/petsc/
> petsc/commits/5675c17
> > DM: Added DMAdaptor
> >
> >   src/snes/utils/dmadapt.c:610
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-linux-c89_el6.log]
> >   /sandbox/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:610:5:
> error: expected expression before '/' token
> >
> >   src/snes/utils/dmadapt.c:611
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-linux-c89_el6.log]
> >   /sandbox/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:611:5:
> error: expected expression before '/' token
> >
> >   src/snes/utils/dmadapt.c:710
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-osx-10.6-cxx-pkgs-opt_ipro.log]
> >   /Users/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:710:16:
> error: no matching function for call to 'PetscFECreateDefault'
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-freebsd-cxx-pkgs-opt_wii.log]
> >   /usr/home/balay/petsc.clone-2/src/snes/utils/dmadapt.c:710:110:
> error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> const char*, PetscInt, _p_PetscFE**)'
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-linux-xsdk-dbg_cg.log]
> >   /sandbox/petsc/petsc.clone/src/snes/utils/dmadapt.c:710:110:
> error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> const char*, PetscInt, _p_PetscFE**)'
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-osx-10.6-cxx-cmplx-pkgs-dbg_ipro.log]
> >   /Users/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:16:
> error: no matching function for call to 'PetscFECreateDefault'
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-linux-opt-cxx-quad_grind.log]
> >   /sandbox/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:110:
> error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> const char*, PetscInt, _p_PetscFE**)'
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-freebsd-cxx-cmplx-pkgs-dbg_wii.log]
> >   /usr/home/balay/petsc.clone/src/snes/utils/dmadapt.c:710:110:
> error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> const char*, PetscInt, _p_PetscFE**)'
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-linux-pkgs-cxx-mlib_el6.log]
> >   /sandbox/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:110:
> error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> const char*, PetscInt, _p_PetscFE**)'
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-linux-cxx-cmplx-pkgs-64idx_churn.log]
> >   /sandbox/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:710:16:
> error: no matching function for call to 'PetscFECreateDefault'
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-linux-xsdk-latest_cg.log]
> >   /sandbox/petsc/petsc.clone-4/src/snes/utils/dmadapt.c:710:110:
> error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> const char*, PetscInt, _p_PetscFE**)'
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 2017/10/17/build_next_arch-linux-gcc-absoft_crank.log]
> >   /sandbox/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:110:
> error: cannot convert 'bool' to 'PetscBool' for argument '4' to
> 'PetscErrorCode PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool,
> const char*, PetscInt, _p_PetscFE**)'
> >
> >   src/snes/utils/dmadapt.c:712
> > [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/
> 

Re: [petsc-dev] [petsc-checkbuilds] PETSc blame digest (next) 2017-10-17

[Satish Balay]

Also - this stuff got into master..

Satish

===

http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/16/examples_full_master.log


not ok diff-ts_tutorials-ex45_2d_q1_r5
#   2,6c2,7
#   < 0 SNES Function norm 2.27683 
#   <   0 KSP Residual norm 6.20756 
#   <   1 KSP Residual norm < 1.e-11
#   < Linear solve converged due to CONVERGED_RTOL iterations 1
#   < 1 SNES Function norm < 1.e-11
#   ---
#   > 0 SNES Function norm 3.20224 
#   >   0 KSP Residual norm 12.6037 
#   >   1 KSP Residual norm 0.000147846 
#   >   2 KSP Residual norm 5.76308e-06 
#   > Linear solve converged due to CONVERGED_RTOL iterations 2
#   > 1 SNES Function norm 9.2212e-06 
#   9,13c10,15
#   < 0 SNES Function norm 2.27316 
#   <   0 KSP Residual norm 6.20174 
#   <   1 KSP Residual norm < 1.e-11
#   < Linear solve converged due to CONVERGED_RTOL iterations 1
#   < 1 SNES Function norm < 1.e-11
#   ---
#   > 0 SNES Function norm 3.20094 
#   >   0 KSP Residual norm 12.6006 
#   >   1 KSP Residual norm 0.00015813 
#   >   2 KSP Residual norm 2.04266e-05 
#   > Linear solve converged due to CONVERGED_RTOL iterations 2
#   > 1 SNES Function norm 1.04029e-05 
#   16,20c18,23
#   < 0 SNES Function norm 2.27313 
#   <   0 KSP Residual norm 6.20054 
#   <   1 KSP Residual norm < 1.e-11
#   < Linear solve converged due to CONVERGED_RTOL iterations 1
#   < 1 SNES Function norm < 1.e-11
#   ---
#   > 0 SNES Function norm 3.20094 
#   >   0 KSP Residual norm 12.6 
#   >   1 KSP Residual norm 0.000144921 
#   >   2 KSP Residual norm 2.52735e-05 
#   > Linear solve converged due to CONVERGED_RTOL iterations 2
#   > 1 SNES Function norm 1.19728e-05 
#   23,27c26,31
#   < 0 SNES Function norm 2.27313 
#   <   0 KSP Residual norm 6.20018 
#   <   1 KSP Residual norm < 1.e-11
#   < Linear solve converged due to CONVERGED_RTOL iterations 1
#   < 1 SNES Function norm < 1.e-11
#   ---
#   > 0 SNES Function norm 3.20094 
#   >   0 KSP Residual norm 12.5999 
#   >   1 KSP Residual norm 0.000160778 
#   >   2 KSP Residual norm 6.78628e-06 
#   > Linear solve converged due to CONVERGED_RTOL iterations 2
#   > 1 SNES Function norm 1.25787e-05 
#   30,34c34,39
#   < 0 SNES Function norm 2.27313 
#   <   0 KSP Residual norm 6.20006 
#   <   1 KSP Residual norm < 1.e-11
#   < Linear solve converged due to CONVERGED_RTOL iterations 1
#   < 1 SNES Function norm < 1.e-11
#   ---
#   > 0 SNES Function norm 3.20094 
#   >   0 KSP Residual norm 12.5998 
#   >   1 KSP Residual norm 0.000171258 
#   >   2 KSP Residual norm 2.25284e-05 
#   > Linear solve converged due to CONVERGED_RTOL iterations 2
#   > 1 SNES Function norm 1.42435e-05 
#   37,41c42,47
#   < 0 SNES Function norm 2.27313 
#   <   0 KSP Residual norm 6.20002 
#   <   1 KSP Residual norm < 1.e-11
#   < Linear solve converged due to CONVERGED_RTOL iterations 1
#   < 1 SNES Function norm < 1.e-11
#   ---
#   > 0 SNES Function norm 3.20094 
#   >   0 KSP Residual norm 12.5998 
#   >   1 KSP Residual norm 0.000169346 
#   >   2 KSP Residual norm 1.85189e-05 
#   > Linear solve converged due to CONVERGED_RTOL iterations 2
#   > 1 SNES Function norm 1.61832e-05 
#   44,48c50,55
#   < 0 SNES Function norm 2.27313 
#   <   0 KSP Residual norm 6.20001 
#   <   1 KSP Residual norm < 1.e-11
#   < Linear solve converged due to CONVERGED_RTOL iterations 1
#   < 1 SNES Function norm < 1.e-11
#   ---
#   > 0 SNES Function norm 3.20094 
#   >   0 KSP Residual norm 12.5998 
#   >   1 KSP Residual norm 0.000172714 
#   >   2 KSP Residual norm 1.6254e-05 
#   > Linear solve converged due to CONVERGED_RTOL iterations 2
#   > 1 SNES Function norm 1.68399e-05 
#   51,55c58,63
#   < 0 SNES Function norm 2.27313 
#   <   0 KSP Residual norm 6.2 
#   <   1 KSP Residual norm < 1.e-11
#   < Linear solve converged due to CONVERGED_RTOL iterations 1
#   < 1 SNES Function norm < 1.e-11
#   ---
#   > 0 SNES Function norm 3.20094 
#   >   0 KSP Residual norm 12.5997 
#   >   1 KSP Residual norm 0.000166623 
#   >   2 KSP Residual norm 1.95142e-05 
#   > Linear solve converged due to CONVERGED_RTOL iterations 2
#   > 1 SNES Function norm 1.80395e-05 
#   58,62c66,71
#   < 0 SNES Function norm 2.27313 
#   <   0 KSP Residual norm 6.2 
#

Re: [petsc-dev] Threadsafe for Matrix assemby in Petsc?

[Mark Adams]

There is no support for threaded matrix assembly in PETSc. Here is a recent
email thread on the issue:

https://mail.google.com/mail/u/0/#search/label%3Apetsc+thread+assembly/15f10d078e9ea8e7

So you pretty much have to deal with race conditions yourself. There are
several failure modes with threads:

1) Off processor entries are stashed in a global data structure for a
scatter/gather stage during matrix finalize. A simple fix for this is to
have every processor compute all elements that touch its vertices
(overlapping element decomposition) and then have PETSc ignore off
processor entries. This is how I do it to simply avoid communication with
redundant computation. It also avoid synchronization.

2) The 1D array data structure is reconstructed when the data spills the
preallocated memory. Solution: allocate memory exactly.

3) Just normal race conditions. Coloring is the basic approach to this
problem, although I think Jed thinks this is not sufficient in PETSc as is.

4) unknown unknowns.

Mark


On Tue, Oct 17, 2017 at 3:24 PM, Yoon, Eisung  wrote:

> Hi Mark,
>
>
>
> Seegyoung here in SCOREC is looking for a way to assemble elements for a
> global matrix with PETSC using threads. I’m aware that you have installed
> thread-safe PETSC into NERSC system for XGC and you know details about
> PETSC.
>
> So could you tell us if PETSC supports thread-safe global matrix assembly
> with thread-safe version of PETSC, and some details, please?
>
> Seegyoung might describe more details about her problem later.
>
>
>
> Thank you.
>
> Best,
>
> ES
>

Re: [petsc-dev] [petsc-checkbuilds] PETSc blame digest (next) 2017-10-17

[Satish Balay]

This has been broken for a couple of days. I pushed a fix [hoping
builds don't abort anymore] - but there are more things to fix..

Satish

On Tue, 17 Oct 2017, PETSc checkBuilds wrote:

> 
> 
> Dear PETSc developer,
> 
> This email contains listings of contributions attributed to you by
> `git blame` that caused compiler errors or warnings in PETSc automated
> testing.  Follow the links to see the full log files. Please attempt to fix
> the issues promptly or let us know at petsc-dev@mcs.anl.gov if you are unable
> to resolve the issues.
> 
> Thanks,
>   The PETSc development team
> 
> 
> 
> warnings attributed to commit 
> https://bitbucket.org/petsc/petsc/commits/5675c17
> DM: Added DMAdaptor
> 
>   src/snes/utils/dmadapt.c:610
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-linux-c89_el6.log]
>   /sandbox/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:610:5: error: 
> expected expression before '/' token
> 
>   src/snes/utils/dmadapt.c:611
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-linux-c89_el6.log]
>   /sandbox/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:611:5: error: 
> expected expression before '/' token
> 
>   src/snes/utils/dmadapt.c:710
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-osx-10.6-cxx-pkgs-opt_ipro.log]
>   /Users/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:710:16: error: no 
> matching function for call to 'PetscFECreateDefault'
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-freebsd-cxx-pkgs-opt_wii.log]
>   /usr/home/balay/petsc.clone-2/src/snes/utils/dmadapt.c:710:110: error: 
> cannot convert 'bool' to 'PetscBool' for argument '4' to 'PetscErrorCode 
> PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool, const char*, 
> PetscInt, _p_PetscFE**)'
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-linux-xsdk-dbg_cg.log]
>   /sandbox/petsc/petsc.clone/src/snes/utils/dmadapt.c:710:110: error: 
> cannot convert 'bool' to 'PetscBool' for argument '4' to 'PetscErrorCode 
> PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool, const char*, 
> PetscInt, _p_PetscFE**)'
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-osx-10.6-cxx-cmplx-pkgs-dbg_ipro.log]
>   /Users/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:16: error: no 
> matching function for call to 'PetscFECreateDefault'
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-linux-opt-cxx-quad_grind.log]
>   /sandbox/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:110: error: 
> cannot convert 'bool' to 'PetscBool' for argument '4' to 'PetscErrorCode 
> PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool, const char*, 
> PetscInt, _p_PetscFE**)'
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-freebsd-cxx-cmplx-pkgs-dbg_wii.log]
>   /usr/home/balay/petsc.clone/src/snes/utils/dmadapt.c:710:110: error: 
> cannot convert 'bool' to 'PetscBool' for argument '4' to 'PetscErrorCode 
> PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool, const char*, 
> PetscInt, _p_PetscFE**)'
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-linux-pkgs-cxx-mlib_el6.log]
>   /sandbox/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:110: error: 
> cannot convert 'bool' to 'PetscBool' for argument '4' to 'PetscErrorCode 
> PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool, const char*, 
> PetscInt, _p_PetscFE**)'
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-linux-cxx-cmplx-pkgs-64idx_churn.log]
>   /sandbox/petsc/petsc.clone-2/src/snes/utils/dmadapt.c:710:16: error: no 
> matching function for call to 'PetscFECreateDefault'
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-linux-xsdk-latest_cg.log]
>   /sandbox/petsc/petsc.clone-4/src/snes/utils/dmadapt.c:710:110: error: 
> cannot convert 'bool' to 'PetscBool' for argument '4' to 'PetscErrorCode 
> PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool, const char*, 
> PetscInt, _p_PetscFE**)'
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-linux-gcc-absoft_crank.log]
>   /sandbox/petsc/petsc.clone-3/src/snes/utils/dmadapt.c:710:110: error: 
> cannot convert 'bool' to 'PetscBool' for argument '4' to 'PetscErrorCode 
> PetscFECreateDefault(DM, PetscInt, PetscInt, PetscBool, const char*, 
> PetscInt, _p_PetscFE**)'
> 
>   src/snes/utils/dmadapt.c:712
> 
> [http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/17/build_next_arch-freebsd-cxx-cmplx-pkgs-dbg_wii.log]
>   /usr/home/balay/petsc.clone/src/snes/utils/dmadapt.c:712:104: error: 
> cannot convert 'bool' to 'PetscBool' for argument '4' to 'PetscErrorCode 
> 

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Barry Smith  writes:

>> On Oct 17, 2017, at 10:41 AM, Jed Brown  wrote:
>> 
>> Barry Smith  writes:
>> 
>> My
>> preference is for TSCreateAdjointTS() which the user can customize to
>> use arbitrary time integration methods (thus becoming a
>> "continuous-in-time" adjoint) and to use different spatial
>> discretizations (for consistency with the adjoint PDE).
> 
>  I understand, I think TSCreateAdjointTS() is not the right model. I 
> think that the correct model is, for example
> 
> PetscSensiCreate()
> PetscSensiSetObjective or whatever etc etc
> PetscSensiGetIntegrator(,)
> PetscSensiGetAdjointIntegrator(,).
 
 How would PetscSensiGetAdjointIntegrator be implemented?  Since it needs
 to be able to create a discrete adjoint (which depends on the forward TS
 implementation), the control flow *must* go through TS.  What does that
 interface look like?
>>> 
>>>   I am not sure what you mean. 
>>> 
>>>   PetscSensi will contain the information about objective function etc as 
>>> well as a pointer to the trajectory history and (of course) a pointer to 
>>> the original TS. All the  auxiliary vectors etc related to sensitivities. 
>>> 
>>>   The adjoint TS that is returned has the calls need to get access to the 
>>> forward solution that it needs etc from the PetscSensi object (which likely 
>>> will call the forward TS for missing timestep values etc. But the calls go 
>>> through the PetscSensi object from the adjoint TS, they don't go directly 
>>> to the forward TS.)
>> 
>> To be concrete with TSRK, PetscSensi does not have access to the Butcher
>> table for the forward method because that only exists in
>> src/ts/impls/explicit/rk/rk.c.  Since the discrete adjoint depends on
>> this detail of the forward method, there MUST be some code somewhere
>> that depends on the forward Butcher table.  Since that data ONLY EXISTS
>> in rk.c, whatever it is that PetscSensi does, there must be some code in
>> rk.c to define the adjoint method.
>
>Of course, that would be incorporated in what is an essentially a 
> TSRKADJOINT implementation that does the adjoint integration. Note that (1) 
> this is simpler than the regular TSRK because it only needs to handle linear 
> with given time steps (2) it also computes along sensitivity information with 
> calls to TSAdjointComputeDRDPFunction() etc
>
>Here is what the code looks like now. Just needs minor refactoring to make 
> it match my proposal.

Right, so this function needs to call the dF/dp (RHSJacobian), dr/dy,
and dr/dp functions.  You're proposing that those would be set on
PetscSensi, but this code below will continue to live in rk.c after
refactoring (rather than in PetscSensi).

Since PetscSensi is not inside TS, there needs to be a public TS
interface ("developer level" if you like, but still public).  I want to
spec out that essential bit of public TS interface.

Note that I proposed putting the sensitivity update stuff below into
post-stage functions.  If that is possible, the bulk of the function
could actually live in PetscSensi.  Or it could be packaged with TS and
reusable across TS implementations.  I still want to start by spec'ing
out what the TS public interface needs to able to support this
PetscSensi object.

> static PetscErrorCode TSAdjointStep_RK(TS ts)
> {
>   TS_RK   *rk   = (TS_RK*)ts->data;
>   RKTableautab  = rk->tableau;
>   const PetscInt   s= tab->s;
>   const PetscReal *A = tab->A,*b = tab->b,*c = tab->c;
>   PetscScalar *w= rk->work;
>   Vec *Y= rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = 
> rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp;
>   PetscInt i,j,nadj;
>   PetscRealt = ts->ptime;
>   PetscErrorCode   ierr;
>   PetscRealh = ts->time_step;
>   Mat  J,Jp;
>
>   PetscFunctionBegin;
>   rk->status = TS_STEP_INCOMPLETE;
>   for (i=s-1; i>=0; i--) {
> rk->stage_time = t + h*(1.0-c[i]);
> for (nadj=0; nadjnumcost; nadj++) {
>   ierr = VecCopy(ts->vecs_sensi[nadj],VecSensiTemp[nadj]);CHKERRQ(ierr);
>   ierr = VecScale(VecSensiTemp[nadj],-h*b[i]);CHKERRQ(ierr);
>   for (j=i+1; j ierr = 
> VecAXPY(VecSensiTemp[nadj],-h*A[j*s+i],VecDeltaLam[nadj*s+j]);CHKERRQ(ierr);
>   }
> }
> /* Stage values of lambda */
> ierr = TSGetRHSJacobian(ts,,,NULL,NULL);CHKERRQ(ierr);
> ierr = TSComputeRHSJacobian(ts,rk->stage_time,Y[i],J,Jp);CHKERRQ(ierr);
> if (ts->vec_costintegral) {
>   ierr = 
> TSAdjointComputeDRDYFunction(ts,rk->stage_time,Y[i],ts->vecs_drdy);CHKERRQ(ierr);
> }
> for (nadj=0; nadjnumcost; nadj++) {
>   ierr = 
> MatMultTranspose(J,VecSensiTemp[nadj],VecDeltaLam[nadj*s+i]);CHKERRQ(ierr);
>   if (ts->vec_costintegral) {
> ierr = 
> VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);CHKERRQ(ierr);
>   }
>   

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Barry Smith]

> On Oct 17, 2017, at 10:41 AM, Jed Brown  wrote:
> 
> Barry Smith  writes:
> 
> My
> preference is for TSCreateAdjointTS() which the user can customize to
> use arbitrary time integration methods (thus becoming a
> "continuous-in-time" adjoint) and to use different spatial
> discretizations (for consistency with the adjoint PDE).
 
  I understand, I think TSCreateAdjointTS() is not the right model. I think 
 that the correct model is, for example
 
 PetscSensiCreate()
 PetscSensiSetObjective or whatever etc etc
 PetscSensiGetIntegrator(,)
 PetscSensiGetAdjointIntegrator(,).
>>> 
>>> How would PetscSensiGetAdjointIntegrator be implemented?  Since it needs
>>> to be able to create a discrete adjoint (which depends on the forward TS
>>> implementation), the control flow *must* go through TS.  What does that
>>> interface look like?
>> 
>>   I am not sure what you mean. 
>> 
>>   PetscSensi will contain the information about objective function etc as 
>> well as a pointer to the trajectory history and (of course) a pointer to the 
>> original TS. All the  auxiliary vectors etc related to sensitivities. 
>> 
>>   The adjoint TS that is returned has the calls need to get access to the 
>> forward solution that it needs etc from the PetscSensi object (which likely 
>> will call the forward TS for missing timestep values etc. But the calls go 
>> through the PetscSensi object from the adjoint TS, they don't go directly to 
>> the forward TS.)
> 
> To be concrete with TSRK, PetscSensi does not have access to the Butcher
> table for the forward method because that only exists in
> src/ts/impls/explicit/rk/rk.c.  Since the discrete adjoint depends on
> this detail of the forward method, there MUST be some code somewhere
> that depends on the forward Butcher table.  Since that data ONLY EXISTS
> in rk.c, whatever it is that PetscSensi does, there must be some code in
> rk.c to define the adjoint method.

   Of course, that would be incorporated in what is an essentially a 
TSRKADJOINT implementation that does the adjoint integration. Note that (1) 
this is simpler than the regular TSRK because it only needs to handle linear 
with given time steps (2) it also computes along sensitivity information with 
calls to TSAdjointComputeDRDPFunction() etc

   Here is what the code looks like now. Just needs minor refactoring to make 
it match my proposal.

static PetscErrorCode TSAdjointStep_RK(TS ts)
{
  TS_RK   *rk   = (TS_RK*)ts->data;
  RKTableautab  = rk->tableau;
  const PetscInt   s= tab->s;
  const PetscReal *A = tab->A,*b = tab->b,*c = tab->c;
  PetscScalar *w= rk->work;
  Vec *Y= rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = 
rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp;
  PetscInt i,j,nadj;
  PetscRealt = ts->ptime;
  PetscErrorCode   ierr;
  PetscRealh = ts->time_step;
  Mat  J,Jp;

  PetscFunctionBegin;
  rk->status = TS_STEP_INCOMPLETE;
  for (i=s-1; i>=0; i--) {
rk->stage_time = t + h*(1.0-c[i]);
for (nadj=0; nadjnumcost; nadj++) {
  ierr = VecCopy(ts->vecs_sensi[nadj],VecSensiTemp[nadj]);CHKERRQ(ierr);
  ierr = VecScale(VecSensiTemp[nadj],-h*b[i]);CHKERRQ(ierr);
  for (j=i+1; jstage_time,Y[i],J,Jp);CHKERRQ(ierr);
if (ts->vec_costintegral) {
  ierr = 
TSAdjointComputeDRDYFunction(ts,rk->stage_time,Y[i],ts->vecs_drdy);CHKERRQ(ierr);
}
for (nadj=0; nadjnumcost; nadj++) {
  ierr = 
MatMultTranspose(J,VecSensiTemp[nadj],VecDeltaLam[nadj*s+i]);CHKERRQ(ierr);
  if (ts->vec_costintegral) {
ierr = 
VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);CHKERRQ(ierr);
  }
}

/* Stage values of mu */
if(ts->vecs_sensip) {
  ierr = 
TSAdjointComputeRHSJacobian(ts,rk->stage_time,Y[i],ts->Jacp);CHKERRQ(ierr);
  if (ts->vec_costintegral) {
ierr = 
TSAdjointComputeDRDPFunction(ts,rk->stage_time,Y[i],ts->vecs_drdp);CHKERRQ(ierr);
  }

  for (nadj=0; nadjnumcost; nadj++) {
ierr = 
MatMultTranspose(ts->Jacp,VecSensiTemp[nadj],VecDeltaMu[nadj*s+i]);CHKERRQ(ierr);
if (ts->vec_costintegral) {
  ierr = 
VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);CHKERRQ(ierr);
}
  }
}
  }

  for (j=0; jvecs_sensi[nadj],s,w,[nadj*s]);CHKERRQ(ierr);
if(ts->vecs_sensip) {
  ierr = 
VecMAXPY(ts->vecs_sensip[nadj],s,w,[nadj*s]);CHKERRQ(ierr);
}
  }
  rk->status = TS_STEP_COMPLETE;
  PetscFunctionReturn(0);
}

> 
> What is that code and what is the call stack to reach it?

Now you ask, 

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Stefano Zampini  writes:

> attached the patch for ts ex2.c to reproduce the problem.
> It runs 30 times the same TSSolve with a TS that has a RHSJacobian and
> BEULER as a solver.

Thanks.

https://bitbucket.org/petsc/petsc/commits/1e3d8eccaccf8d08bd0d0444380c7e269e8f696e

> This are the times for each solve on my workstation
>
> [szampini@localhost tutorials]$ ./ex2
> 0.164751
> 0.1406
> 0.148703
> 0.158375
> 0.169456
> 0.172364
> 0.179924
> 0.186967
> 0.195149
> 0.202559
> 0.214942
> 0.217264
> 0.225801
> 0.232413
> 0.240205
> 0.249077
> 0.256062
> 0.262738
> 0.270429
> 0.277631
> 0.285879
> 0.292678
> 0.299828
> 0.307717
> 0.315533
> 0.323094
> 0.330388
> 0.338611
> 0.345624
> 0.35268
>
> You can see that timings (for the same exact solve) increase. If I comment
> out DMCoarsenHookAdd in SNESSetUpMatrices
>
> [szampini@localhost tutorials]$ git diff
> /home/szampini/src/petsc/src/snes/interface/
> diff --git a/src/snes/interface/snes.c b/src/snes/interface/snes.c
> index d366a2cc69..cdb7571869 100644
> --- a/src/snes/interface/snes.c
> +++ b/src/snes/interface/snes.c
> @@ -642,7 +642,7 @@ PetscErrorCode SNESSetUpMatrices(SNES snes)
>  KSP ksp;
>  ierr = SNESGetKSP(snes,);CHKERRQ(ierr);
>  ierr =
> KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);CHKERRQ(ierr);
> -ierr =
> DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);CHKERRQ(ierr);
> +//ierr =
> DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);CHKERRQ(ierr);
>}
>PetscFunctionReturn(0);
>  }
>
> this is what I get, constant timings as expected
>
> [szampini@localhost tutorials]$ ./ex2
> 0.154417
> 0.125632
> 0.125743
> 0.126303
> 0.125648
> 0.125593
> 0.125647
> 0.125656
> 0.125705
> 0.125705
> 0.125737
> 0.125648
> 0.125592
> 0.125568
> 0.125578
> 0.125678
> 0.125628
> 0.126056
> 0.125725
> 0.127705
> 0.126028
> 0.125671
> 0.125675
> 0.125662
> 0.125688
> 0.125767
> 0.12574
> 0.125604
> 0.125568
> 0.1256

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Lawrence Mitchell]

> On 17 Oct 2017, at 17:35, Jed Brown  wrote:
> 
>> 
>> So my initial concern was that if I do:
>> 
>> SNESSetDM(snes, dm);
>> SNESSolve(snes);
>> SNESDestroy();
>> 
>> Then dm still has stale pointers to snes (in the context of the
>> restrict hook).  
> 
> Thanks, this is fixable.  Can someone make a test case for me?

I'll try and cook one up tomorrow.

Cheers,

Lawrence

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Lawrence Mitchell  writes:

>> On 17 Oct 2017, at 15:47, Jed Brown  wrote:
>> 
>>> c.f. also the petsc-maint discussion last October in the thread:
>>> 
>>> "Segfault when DM is reused in two SNESes with multigrid"
>> 
>> What should I do?  That thread petered out when nobody could suggest a
>> criteria to distinguish reuse of a DM in related SNES objects (e.g.,
>> nonlinear preconditioning or nonlinear multigrid) from reuse with intent
>> to solve different problems.
>
> So my initial concern was that if I do:
>
> SNESSetDM(snes, dm);
> SNESSolve(snes);
> SNESDestroy();
>
> Then dm still has stale pointers to snes (in the context of the
> restrict hook).  

Thanks, this is fixable.  Can someone make a test case for me?

> These aren't cleaned up until DMDestroy.  It was surprising to me that
> SNESDestroy did not clean up all the things that SNESCreate/Solve set
> up.
>
> I don't care about the case of:
>
> SNESSetDM(snes1, dm);
> SNESSetDM(snes2, dm);
>
> I work around this issue by making new DMs, so it's OK if nothing changes in 
> PETSc.  
>
> Lawrence

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Stefano Zampini]

attached the patch for ts ex2.c to reproduce the problem.
It runs 30 times the same TSSolve with a TS that has a RHSJacobian and
BEULER as a solver.
This are the times for each solve on my workstation

[szampini@localhost tutorials]$ ./ex2
0.164751
0.1406
0.148703
0.158375
0.169456
0.172364
0.179924
0.186967
0.195149
0.202559
0.214942
0.217264
0.225801
0.232413
0.240205
0.249077
0.256062
0.262738
0.270429
0.277631
0.285879
0.292678
0.299828
0.307717
0.315533
0.323094
0.330388
0.338611
0.345624
0.35268

You can see that timings (for the same exact solve) increase. If I comment
out DMCoarsenHookAdd in SNESSetUpMatrices

[szampini@localhost tutorials]$ git diff
/home/szampini/src/petsc/src/snes/interface/
diff --git a/src/snes/interface/snes.c b/src/snes/interface/snes.c
index d366a2cc69..cdb7571869 100644
--- a/src/snes/interface/snes.c
+++ b/src/snes/interface/snes.c
@@ -642,7 +642,7 @@ PetscErrorCode SNESSetUpMatrices(SNES snes)
 KSP ksp;
 ierr = SNESGetKSP(snes,);CHKERRQ(ierr);
 ierr =
KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);CHKERRQ(ierr);
-ierr =
DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);CHKERRQ(ierr);
+//ierr =
DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }

this is what I get, constant timings as expected

[szampini@localhost tutorials]$ ./ex2
0.154417
0.125632
0.125743
0.126303
0.125648
0.125593
0.125647
0.125656
0.125705
0.125705
0.125737
0.125648
0.125592
0.125568
0.125578
0.125678
0.125628
0.126056
0.125725
0.127705
0.126028
0.125671
0.125675
0.125662
0.125688
0.125767
0.12574
0.125604
0.125568
0.1256


patch
Description: Binary data

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Lawrence Mitchell]

> On 17 Oct 2017, at 15:47, Jed Brown  wrote:
> 
>> c.f. also the petsc-maint discussion last October in the thread:
>> 
>> "Segfault when DM is reused in two SNESes with multigrid"
> 
> What should I do?  That thread petered out when nobody could suggest a
> criteria to distinguish reuse of a DM in related SNES objects (e.g.,
> nonlinear preconditioning or nonlinear multigrid) from reuse with intent
> to solve different problems.

So my initial concern was that if I do:

SNESSetDM(snes, dm);
SNESSolve(snes);
SNESDestroy();

Then dm still has stale pointers to snes (in the context of the restrict hook). 
 These aren't cleaned up until DMDestroy.  It was surprising to me that 
SNESDestroy did not clean up all the things that SNESCreate/Solve set up.

I don't care about the case of:

SNESSetDM(snes1, dm);
SNESSetDM(snes2, dm);

I work around this issue by making new DMs, so it's OK if nothing changes in 
PETSc.  

Lawrence

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Barry Smith  writes:

>   Ok, I think it plus the DMTS are problematically because, frankly,
>   only you understand them and anything in a package that has subtle
>   complexities in it that only one person understands is really bad
>   (you kind of admit this by saying you are the only one who can
>   really do a new refactorization).

I don't think it's much different from other parts of PETSc.  How many
people could jump in and refactor PCGAMG or PCFIELDSPLIT or SNESNASM or
VecScatter or DMPlex?  I think it would take any of us a few hours (or
more) of reading, understanding each use case, and sketching how a
replacement can provide each of those use cases.  This is vulnerable to
Second System Syndrome and someone to spend time making a several
thousand line patch with sweeping changes that wasn't carefully thought
out (and preferably incremental), nor would I want to review such a
beast.

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Stefano Zampini  writes:

>> >> have you thought about how to fix this?
>> >>
>> >> https://lists.mcs.anl.gov/pipermail/petsc-dev/2017-October/021387.html
>>
>> I feel like I'm missing context.  I know I failed to follow up in a
>> thread at some time in the past, but I don't really understand the
>> problem or what is being attempted.  What is the failing test case?
>>
>>
> There's no failing test; however, you can see considerable slow downs when
> calling TSSolve on the same ts multiple times. Just copying and paste the
> original mail message below: it seems clear to me the code path. However,
> what it's not clear is why we need to call DMCoarsenHookAdd so many times.
> Cannot be set once?
>
> Every time you call TSGetRHSMats_Private, TSGetIJacobian gets called,
> which in turns calls SNESSetUpMatrices, that adds the entry in the linked
> list of hooks via DMCoarsenHookAdd. This causes a considerable slow down
> when using the RHS interface with an implicit solver, as you can see from
> running the patched code.

What is "the patched code"?  I can't test whether a fix works without
being able to run something.  Since you have encountered this in actual
code, I'm asking you to provide a test case that I can run to see this
effect and to confirm if I fix it.  Put it in a branch and make a PR
with me responsible.  It's work you've clearly done in some form and I
shouldn't have to reproduce, especially when I don't get huge blocks of
time to work on PETSc and have plenty of research-critical tasks (in
PETSc and otherwise) that I'm unable to get to.

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Stefano Zampini]

> >> have you thought about how to fix this?
> >>
> >> https://lists.mcs.anl.gov/pipermail/petsc-dev/2017-October/021387.html
>
> I feel like I'm missing context.  I know I failed to follow up in a
> thread at some time in the past, but I don't really understand the
> problem or what is being attempted.  What is the failing test case?
>
>
There's no failing test; however, you can see considerable slow downs when
calling TSSolve on the same ts multiple times. Just copying and paste the
original mail message below: it seems clear to me the code path. However,
what it's not clear is why we need to call DMCoarsenHookAdd so many times.
Cannot be set once?

Every time you call TSGetRHSMats_Private, TSGetIJacobian gets called,
which in turns calls SNESSetUpMatrices, that adds the entry in the linked
list of hooks via DMCoarsenHookAdd. This causes a considerable slow down
when using the RHS interface with an implicit solver, as you can see from
running the patched code.

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Barry Smith  writes:

>> On Oct 17, 2017, at 10:00 AM, Jed Brown  wrote:
>> 
>> Barry Smith  writes:
>> 
 On Oct 17, 2017, at 9:20 AM, Jed Brown  wrote:
 
 Barry Smith  writes:
 
>> On Oct 17, 2017, at 7:43 AM, Stefano Zampini  
>> wrote:
>> 
>> 
>> 
>> 
>>  What happens if you call TSCreateAdjointsTS() on a TS obtained with 
>> TSCreateAdjointsTS()? Is the resulting TS useful?
>> 
>> 
>> I don't think so.
> 
>  Hmm, this is rather bothersome. So the new TS that comes out is not
>  really a full-featured TS (using your new extended TS object that
>  knows about adjoints)? 
 
 Is it useful to precondition Newton with Newton preconditioned with
 Newton?  
 
>  It is like a subclass of TS that can only do timestepping? This
>  goes back to my concern that TS shouldn't be overloaded with all
>  the adjoint sensitivity stuff, there should be a higher level
>  PetscSensi object that does sensitivities and it contains multiple
>  TS for the ODE/DAE integrator and one for "reverse" mode.
 
 Let's not mix together Stefano's TSComputeObjectiveAndGradient (which
 could easily go into a helper PetscSensi) with the control flow needed
 to create/integrate an adjoint system.
>>> 
>>>   Jed,
>>> 
>>>I am not mixing them!  I am not even considering the 
>>> TSComputeObjectiveAndGradient() at all I am only considering
>>> the "the control flow needed to create/integrate an adjoint system." My 
>>> PetscSensi is not for TSComputeObjectiveAndGradient level stuff.
>> 
>> Okay, interesting.
>> 
>>> 
 In my opinion, it must be
 possible to have a comparable/same interface to create a discrete
 adjoint (which depends on the forward implementation) as to create a
 continuous adjoint (with customizable space and time discretization).
>>> 
>>>   Absolutely and fundamentally. No argument at all about this.
>>> 
 
 I think we should postpone discussion of higher level convenience
 functions like TSComputeObjectiveAndGradient or some new organization
 thereof until we agree about the basic representation of adjoints.
>>> 
>>>   Totally agree and that is what I have done. 
>>> 
 My
 preference is for TSCreateAdjointTS() which the user can customize to
 use arbitrary time integration methods (thus becoming a
 "continuous-in-time" adjoint) and to use different spatial
 discretizations (for consistency with the adjoint PDE).
>>> 
>>>   I understand, I think TSCreateAdjointTS() is not the right model. I think 
>>> that the correct model is, for example
>>> 
>>> PetscSensiCreate()
>>> PetscSensiSetObjective or whatever etc etc
>>> PetscSensiGetIntegrator(,)
>>> PetscSensiGetAdjointIntegrator(,).
>> 
>> How would PetscSensiGetAdjointIntegrator be implemented?  Since it needs
>> to be able to create a discrete adjoint (which depends on the forward TS
>> implementation), the control flow *must* go through TS.  What does that
>> interface look like?
>
>I am not sure what you mean. 
>
>PetscSensi will contain the information about objective function etc as 
> well as a pointer to the trajectory history and (of course) a pointer to the 
> original TS. All the  auxiliary vectors etc related to sensitivities. 
>
>The adjoint TS that is returned has the calls need to get access to the 
> forward solution that it needs etc from the PetscSensi object (which likely 
> will call the forward TS for missing timestep values etc. But the calls go 
> through the PetscSensi object from the adjoint TS, they don't go directly to 
> the forward TS.)

To be concrete with TSRK, PetscSensi does not have access to the Butcher
table for the forward method because that only exists in
src/ts/impls/explicit/rk/rk.c.  Since the discrete adjoint depends on
this detail of the forward method, there MUST be some code somewhere
that depends on the forward Butcher table.  Since that data ONLY EXISTS
in rk.c, whatever it is that PetscSensi does, there must be some code in
rk.c to define the adjoint method.

What is that code and what is the call stack to reach it?

>Both proposed models shove one heck of a lot of stuff into TS; I think it 
> would be a better design and better user API if we don't do that; we reserve 
> the TS for being able to integrate something (and that is about all). Maybe 
> it is not possible but I don't see why it is not possible.
>
>Since I don't know the details I asked for the people who do know the 
> details (Hong and Stefano) to try to develop such an API (or clearly 
> articular why it is impossible.)

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Barry Smith]

> On Oct 17, 2017, at 10:30 AM, Jed Brown  wrote:
> 
> Barry Smith  writes:
> 
>>> On Oct 17, 2017, at 10:06 AM, Jed Brown  wrote:
>>> 
>>> Barry Smith  writes:
>>> 
> On Oct 17, 2017, at 9:47 AM, Jed Brown  wrote:
> 
> Lawrence Mitchell  writes:
> 
> 
> When I suggested as a young child that DM be essentially just a function
> space and create a new object for resolution-independent specification
> of a problem (residual and Jacobian functions and related components),
> Barry wanted it to be part of DM to avoid having a new object.  So it's
> part of DM -- make a new DM if you're solving a different problem.
 
 Of course, everything in PETSc is subject to refactorization and it
 may be time to do this refactorization; especially if it can
 dramatically decrease the ugly subtle complexities of the TSDM,
 SNESDM  management. One more public object per solver level is
 probably better than the complexity we have now I do admit.
>>> 
>>> I'm not opposed to refactoring (though it would take me significant time
>>> without distractions),
>> 
>>   Perhaps you don't need to do it all yourself? And we could do it in
>>   multiple stages, like first add the new public objects get them
>>   working with the DM and then later remove the private objects that
>>   exist now?
> 
> Whomever does it will need to spend some time learning/reminding
> themselves of all the functional requirements.  It could be done
> somewhat incrementally, but will need careful reviewing.
> 
>>> but this sort of change would have a lot more
>>> consequences now because we have lots of code depending on it.  
>> 
>>   True
>> 
>>> Is there
>>> a functional reason to refactor now?
>> 
>>Is there ever? As always it is priorities.
> 
> Nonlinear preconditioning and multilevel algorithms were the impetus for
> DMSNES.  Prior to that, we had some kludgy void* hiding a dependency
> loop (DM depending on SNES) and that was definitely not maintainable or
> reasonable to extend to support the new functionality.
> 
> I don't see the current structure as a significant liability and
> wouldn't consider it high priority at this time.

  Ok, I think it plus the DMTS are problematically because, frankly, only you 
understand them and anything in a package that has subtle complexities in it 
that only one person understands is really bad (you kind of admit this by 
saying you are the only one who can really do a new refactorization).

  Barry

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Barry Smith  writes:

>> On Oct 17, 2017, at 10:06 AM, Jed Brown  wrote:
>> 
>> Barry Smith  writes:
>> 
 On Oct 17, 2017, at 9:47 AM, Jed Brown  wrote:
 
 Lawrence Mitchell  writes:
 
 
 When I suggested as a young child that DM be essentially just a function
 space and create a new object for resolution-independent specification
 of a problem (residual and Jacobian functions and related components),
 Barry wanted it to be part of DM to avoid having a new object.  So it's
 part of DM -- make a new DM if you're solving a different problem.
>>> 
>>>  Of course, everything in PETSc is subject to refactorization and it
>>>  may be time to do this refactorization; especially if it can
>>>  dramatically decrease the ugly subtle complexities of the TSDM,
>>>  SNESDM  management. One more public object per solver level is
>>>  probably better than the complexity we have now I do admit.
>> 
>> I'm not opposed to refactoring (though it would take me significant time
>> without distractions),
>
>Perhaps you don't need to do it all yourself? And we could do it in
>multiple stages, like first add the new public objects get them
>working with the DM and then later remove the private objects that
>exist now?

Whomever does it will need to spend some time learning/reminding
themselves of all the functional requirements.  It could be done
somewhat incrementally, but will need careful reviewing.

>> but this sort of change would have a lot more
>> consequences now because we have lots of code depending on it.  
>
>True
>
>> Is there
>> a functional reason to refactor now?
>
> Is there ever? As always it is priorities.

Nonlinear preconditioning and multilevel algorithms were the impetus for
DMSNES.  Prior to that, we had some kludgy void* hiding a dependency
loop (DM depending on SNES) and that was definitely not maintainable or
reasonable to extend to support the new functionality.

I don't see the current structure as a significant liability and
wouldn't consider it high priority at this time.

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Barry Smith]

> On Oct 17, 2017, at 10:06 AM, Jed Brown  wrote:
> 
> Barry Smith  writes:
> 
>>> On Oct 17, 2017, at 9:47 AM, Jed Brown  wrote:
>>> 
>>> Lawrence Mitchell  writes:
>>> 
>>> 
>>> When I suggested as a young child that DM be essentially just a function
>>> space and create a new object for resolution-independent specification
>>> of a problem (residual and Jacobian functions and related components),
>>> Barry wanted it to be part of DM to avoid having a new object.  So it's
>>> part of DM -- make a new DM if you're solving a different problem.
>> 
>>  Of course, everything in PETSc is subject to refactorization and it
>>  may be time to do this refactorization; especially if it can
>>  dramatically decrease the ugly subtle complexities of the TSDM,
>>  SNESDM  management. One more public object per solver level is
>>  probably better than the complexity we have now I do admit.
> 
> I'm not opposed to refactoring (though it would take me significant time
> without distractions),

   Perhaps you don't need to do it all yourself? And we could do it in multiple 
stages, like first add the new public objects get them working with the DM and 
then later remove the private objects that exist now?

> but this sort of change would have a lot more
> consequences now because we have lots of code depending on it.  

   True

> Is there
> a functional reason to refactor now?

Is there ever? As always it is priorities.

  Barry

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Barry Smith]

> On Oct 17, 2017, at 10:00 AM, Jed Brown  wrote:
> 
> Barry Smith  writes:
> 
>>> On Oct 17, 2017, at 9:20 AM, Jed Brown  wrote:
>>> 
>>> Barry Smith  writes:
>>> 
> On Oct 17, 2017, at 7:43 AM, Stefano Zampini  
> wrote:
> 
> 
> 
> 
>  What happens if you call TSCreateAdjointsTS() on a TS obtained with 
> TSCreateAdjointsTS()? Is the resulting TS useful?
> 
> 
> I don't think so.
 
  Hmm, this is rather bothersome. So the new TS that comes out is not
  really a full-featured TS (using your new extended TS object that
  knows about adjoints)? 
>>> 
>>> Is it useful to precondition Newton with Newton preconditioned with
>>> Newton?  
>>> 
  It is like a subclass of TS that can only do timestepping? This
  goes back to my concern that TS shouldn't be overloaded with all
  the adjoint sensitivity stuff, there should be a higher level
  PetscSensi object that does sensitivities and it contains multiple
  TS for the ODE/DAE integrator and one for "reverse" mode.
>>> 
>>> Let's not mix together Stefano's TSComputeObjectiveAndGradient (which
>>> could easily go into a helper PetscSensi) with the control flow needed
>>> to create/integrate an adjoint system.
>> 
>>   Jed,
>> 
>>I am not mixing them!  I am not even considering the 
>> TSComputeObjectiveAndGradient() at all I am only considering
>> the "the control flow needed to create/integrate an adjoint system." My 
>> PetscSensi is not for TSComputeObjectiveAndGradient level stuff.
> 
> Okay, interesting.
> 
>> 
>>> In my opinion, it must be
>>> possible to have a comparable/same interface to create a discrete
>>> adjoint (which depends on the forward implementation) as to create a
>>> continuous adjoint (with customizable space and time discretization).
>> 
>>   Absolutely and fundamentally. No argument at all about this.
>> 
>>> 
>>> I think we should postpone discussion of higher level convenience
>>> functions like TSComputeObjectiveAndGradient or some new organization
>>> thereof until we agree about the basic representation of adjoints.
>> 
>>   Totally agree and that is what I have done. 
>> 
>>> My
>>> preference is for TSCreateAdjointTS() which the user can customize to
>>> use arbitrary time integration methods (thus becoming a
>>> "continuous-in-time" adjoint) and to use different spatial
>>> discretizations (for consistency with the adjoint PDE).
>> 
>>   I understand, I think TSCreateAdjointTS() is not the right model. I think 
>> that the correct model is, for example
>> 
>> PetscSensiCreate()
>> PetscSensiSetObjective or whatever etc etc
>> PetscSensiGetIntegrator(,)
>> PetscSensiGetAdjointIntegrator(,).
> 
> How would PetscSensiGetAdjointIntegrator be implemented?  Since it needs
> to be able to create a discrete adjoint (which depends on the forward TS
> implementation), the control flow *must* go through TS.  What does that
> interface look like?

   I am not sure what you mean. 

   PetscSensi will contain the information about objective function etc as well 
as a pointer to the trajectory history and (of course) a pointer to the 
original TS. All the  auxiliary vectors etc related to sensitivities. 

   The adjoint TS that is returned has the calls need to get access to the 
forward solution that it needs etc from the PetscSensi object (which likely 
will call the forward TS for missing timestep values etc. But the calls go 
through the PetscSensi object from the adjoint TS, they don't go directly to 
the forward TS.)

   Both proposed models shove one heck of a lot of stuff into TS; I think it 
would be a better design and better user API if we don't do that; we reserve 
the TS for being able to integrate something (and that is about all). Maybe it 
is not possible but I don't see why it is not possible.

   Since I don't know the details I asked for the people who do know the 
details (Hong and Stefano) to try to develop such an API (or clearly articular 
why it is impossible.)


  Barry



> 
>> ...
>> 
>>  Note I am not married to the current model of TSSolve/TSAdjointSolve().
>> 
>>  Barry

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Barry Smith  writes:

>> On Oct 17, 2017, at 9:47 AM, Jed Brown  wrote:
>> 
>> Lawrence Mitchell  writes:
>> 
>> 
>> When I suggested as a young child that DM be essentially just a function
>> space and create a new object for resolution-independent specification
>> of a problem (residual and Jacobian functions and related components),
>> Barry wanted it to be part of DM to avoid having a new object.  So it's
>> part of DM -- make a new DM if you're solving a different problem.
>
>   Of course, everything in PETSc is subject to refactorization and it
>   may be time to do this refactorization; especially if it can
>   dramatically decrease the ugly subtle complexities of the TSDM,
>   SNESDM  management. One more public object per solver level is
>   probably better than the complexity we have now I do admit.

I'm not opposed to refactoring (though it would take me significant time
without distractions), but this sort of change would have a lot more
consequences now because we have lots of code depending on it.  Is there
a functional reason to refactor now?

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Barry Smith  writes:

>> On Oct 17, 2017, at 9:20 AM, Jed Brown  wrote:
>> 
>> Barry Smith  writes:
>> 
 On Oct 17, 2017, at 7:43 AM, Stefano Zampini  
 wrote:
 
 
 
 
   What happens if you call TSCreateAdjointsTS() on a TS obtained with 
 TSCreateAdjointsTS()? Is the resulting TS useful?
 
 
 I don't think so.
>>> 
>>>   Hmm, this is rather bothersome. So the new TS that comes out is not
>>>   really a full-featured TS (using your new extended TS object that
>>>   knows about adjoints)? 
>> 
>> Is it useful to precondition Newton with Newton preconditioned with
>> Newton?  
>> 
>>>   It is like a subclass of TS that can only do timestepping? This
>>>   goes back to my concern that TS shouldn't be overloaded with all
>>>   the adjoint sensitivity stuff, there should be a higher level
>>>   PetscSensi object that does sensitivities and it contains multiple
>>>   TS for the ODE/DAE integrator and one for "reverse" mode.
>> 
>> Let's not mix together Stefano's TSComputeObjectiveAndGradient (which
>> could easily go into a helper PetscSensi) with the control flow needed
>> to create/integrate an adjoint system.
>
>Jed,
>
> I am not mixing them!  I am not even considering the 
> TSComputeObjectiveAndGradient() at all I am only considering
> the "the control flow needed to create/integrate an adjoint system." My 
> PetscSensi is not for TSComputeObjectiveAndGradient level stuff.

Okay, interesting.

>
>>  In my opinion, it must be
>> possible to have a comparable/same interface to create a discrete
>> adjoint (which depends on the forward implementation) as to create a
>> continuous adjoint (with customizable space and time discretization).
>
>Absolutely and fundamentally. No argument at all about this.
>
>> 
>> I think we should postpone discussion of higher level convenience
>> functions like TSComputeObjectiveAndGradient or some new organization
>> thereof until we agree about the basic representation of adjoints.
>
>Totally agree and that is what I have done. 
>
>>  My
>> preference is for TSCreateAdjointTS() which the user can customize to
>> use arbitrary time integration methods (thus becoming a
>> "continuous-in-time" adjoint) and to use different spatial
>> discretizations (for consistency with the adjoint PDE).
>
>I understand, I think TSCreateAdjointTS() is not the right model. I think 
> that the correct model is, for example
>
> PetscSensiCreate()
> PetscSensiSetObjective or whatever etc etc
> PetscSensiGetIntegrator(,)
> PetscSensiGetAdjointIntegrator(,).

How would PetscSensiGetAdjointIntegrator be implemented?  Since it needs
to be able to create a discrete adjoint (which depends on the forward TS
implementation), the control flow *must* go through TS.  What does that
interface look like?

> ...
>
>   Note I am not married to the current model of TSSolve/TSAdjointSolve().
>
>   Barry

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Barry Smith]

> On Oct 17, 2017, at 9:47 AM, Jed Brown  wrote:
> 
> Lawrence Mitchell  writes:
> 
> 
> When I suggested as a young child that DM be essentially just a function
> space and create a new object for resolution-independent specification
> of a problem (residual and Jacobian functions and related components),
> Barry wanted it to be part of DM to avoid having a new object.  So it's
> part of DM -- make a new DM if you're solving a different problem.

  Of course, everything in PETSc is subject to refactorization and it may be 
time to do this refactorization; especially if it can dramatically decrease the 
ugly subtle complexities of the TSDM, SNESDM  management. One more public 
object per solver level is probably better than the complexity we have now I do 
admit.

  Barry

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Barry Smith]

> On Oct 17, 2017, at 9:20 AM, Jed Brown  wrote:
> 
> Barry Smith  writes:
> 
>>> On Oct 17, 2017, at 7:43 AM, Stefano Zampini  
>>> wrote:
>>> 
>>> 
>>> 
>>> 
>>>   What happens if you call TSCreateAdjointsTS() on a TS obtained with 
>>> TSCreateAdjointsTS()? Is the resulting TS useful?
>>> 
>>> 
>>> I don't think so.
>> 
>>   Hmm, this is rather bothersome. So the new TS that comes out is not
>>   really a full-featured TS (using your new extended TS object that
>>   knows about adjoints)? 
> 
> Is it useful to precondition Newton with Newton preconditioned with
> Newton?  
> 
>>   It is like a subclass of TS that can only do timestepping? This
>>   goes back to my concern that TS shouldn't be overloaded with all
>>   the adjoint sensitivity stuff, there should be a higher level
>>   PetscSensi object that does sensitivities and it contains multiple
>>   TS for the ODE/DAE integrator and one for "reverse" mode.
> 
> Let's not mix together Stefano's TSComputeObjectiveAndGradient (which
> could easily go into a helper PetscSensi) with the control flow needed
> to create/integrate an adjoint system.

   Jed,

I am not mixing them!  I am not even considering the 
TSComputeObjectiveAndGradient() at all I am only considering
the "the control flow needed to create/integrate an adjoint system." My 
PetscSensi is not for TSComputeObjectiveAndGradient level stuff.


>  In my opinion, it must be
> possible to have a comparable/same interface to create a discrete
> adjoint (which depends on the forward implementation) as to create a
> continuous adjoint (with customizable space and time discretization).

   Absolutely and fundamentally. No argument at all about this.

> 
> I think we should postpone discussion of higher level convenience
> functions like TSComputeObjectiveAndGradient or some new organization
> thereof until we agree about the basic representation of adjoints.

   Totally agree and that is what I have done. 

>  My
> preference is for TSCreateAdjointTS() which the user can customize to
> use arbitrary time integration methods (thus becoming a
> "continuous-in-time" adjoint) and to use different spatial
> discretizations (for consistency with the adjoint PDE).

   I understand, I think TSCreateAdjointTS() is not the right model. I think 
that the correct model is, for example

PetscSensiCreate()
PetscSensiSetObjective or whatever etc etc
PetscSensiGetIntegrator(,)
PetscSensiGetAdjointIntegrator(,).
...

  Note I am not married to the current model of TSSolve/TSAdjointSolve().

  Barry

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Lawrence Mitchell  writes:

>> On 16 Oct 2017, at 09:35, Stefano Zampini  wrote:
>> 
>> Jed,
>> 
>> have you thought about how to fix this?
>> 
>> https://lists.mcs.anl.gov/pipermail/petsc-dev/2017-October/021387.html

I feel like I'm missing context.  I know I failed to follow up in a
thread at some time in the past, but I don't really understand the
problem or what is being attempted.  What is the failing test case?

> c.f. also the petsc-maint discussion last October in the thread:
>
> "Segfault when DM is reused in two SNESes with multigrid"

What should I do?  That thread petered out when nobody could suggest a
criteria to distinguish reuse of a DM in related SNES objects (e.g.,
nonlinear preconditioning or nonlinear multigrid) from reuse with intent
to solve different problems.

When I suggested as a young child that DM be essentially just a function
space and create a new object for resolution-independent specification
of a problem (residual and Jacobian functions and related components),
Barry wanted it to be part of DM to avoid having a new object.  So it's
part of DM -- make a new DM if you're solving a different problem.

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Barry Smith  writes:

>> On Oct 17, 2017, at 7:43 AM, Stefano Zampini  
>> wrote:
>> 
>> 
>> 
>> 
>>What happens if you call TSCreateAdjointsTS() on a TS obtained with 
>> TSCreateAdjointsTS()? Is the resulting TS useful?
>> 
>> 
>> I don't think so.
>
>Hmm, this is rather bothersome. So the new TS that comes out is not
>really a full-featured TS (using your new extended TS object that
>knows about adjoints)? 

Is it useful to precondition Newton with Newton preconditioned with
Newton?  

>It is like a subclass of TS that can only do timestepping? This
>goes back to my concern that TS shouldn't be overloaded with all
>the adjoint sensitivity stuff, there should be a higher level
>PetscSensi object that does sensitivities and it contains multiple
>TS for the ODE/DAE integrator and one for "reverse" mode.

Let's not mix together Stefano's TSComputeObjectiveAndGradient (which
could easily go into a helper PetscSensi) with the control flow needed
to create/integrate an adjoint system.  In my opinion, it must be
possible to have a comparable/same interface to create a discrete
adjoint (which depends on the forward implementation) as to create a
continuous adjoint (with customizable space and time discretization).

I think we should postpone discussion of higher level convenience
functions like TSComputeObjectiveAndGradient or some new organization
thereof until we agree about the basic representation of adjoints.  My
preference is for TSCreateAdjointTS() which the user can customize to
use arbitrary time integration methods (thus becoming a
"continuous-in-time" adjoint) and to use different spatial
discretizations (for consistency with the adjoint PDE).

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Barry Smith]

> On Oct 17, 2017, at 8:56 AM, Stefano Zampini  
> wrote:
> 
> 
> 
> 
> 
>Yes, we are waiting for Hong to provide his proposed API.
> 
>BTW: Stefano has not really submitted a proposed API as I requested, just 
> some minor explanation of his current API and minor changes.
> 
> 
> My proposal would be based on what we are currently discussing, and it will 
> build on the existing one.

  I didn't ask for that. See my previous email. I want one based on a 
PetscSensi object

  Barry

> 
> 
> 
> -- 
> Stefano

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Barry Smith]

> On Oct 17, 2017, at 7:43 AM, Stefano Zampini  
> wrote:
> 
> 
> 
> 
>What happens if you call TSCreateAdjointsTS() on a TS obtained with 
> TSCreateAdjointsTS()? Is the resulting TS useful?
> 
> 
> I don't think so.

   Hmm, this is rather bothersome. So the new TS that comes out is not really a 
full-featured TS (using your new extended TS object that knows about adjoints)? 
It is like a subclass of TS that can only do timestepping? This goes back to my 
concern that TS shouldn't be overloaded with all the adjoint sensitivity stuff, 
there should be a higher level PetscSensi object that does sensitivities and it 
contains multiple TS for the ODE/DAE integrator and one for "reverse" mode.

  Barry

In this model Hong's TSAdjoint... stuff that deals with integration details 
would be in the second sub TS while the stuff that sets up the adjoint problem 
etc would be moved to the PetscSensi object.




> Theoretically, it should be a tangent linear model with a different forcing 
> function
> 
> 
>  
> Barry
> 
> 
> > On Oct 17, 2017, at 12:37 AM, Stefano Zampini  
> > wrote:
> >
> >
> >
> >
> > >
> > >
> > >> In case of multiple objectives, there may be a performance reason to
> > >> amortize evaluation of several at once, though the list interface is
> > >> convenient.  Consider common objectives being quantities like lift and
> > >> drag on different surfaces of a fluids simulation or stress/strain at
> > >> certain critical joints in a structure.  Although these have some
> > >> locality, it's reasonable to assume that state dependence will have
> > >> quickly become global, thus make no attempt to handle sparse
> > >> representations of the adjoint vectors lambda.
> > >>
> > >>
> > > I don't get this comment. Is it related with multi-objective optimization
> > > (e.g. Pareto)?
> >
> > Adjoints are usually preferred any time you are differentiating a small
> > number of output variables with respect to a large number of inputs.  It
> > could be for multi-objective optimization, but it's every bit as
> > relevant for physical sensitivities.
> >
> >
> > If we are not talking about Pareto optimization, and thus we don't need a 
> > separate output from each function, then users can pass a single function 
> > that computes all the quantities they need at the same time. Anyway, I 
> > don't mind having a single callback for multiple functions.
> >
> > >> How are parameters accessed in TSComputeRHSFunction?  It looks like
> > >> they're coming out of the context.  Why should this be different?  (If
> > >> parameters need to go into a Vec, we could do that, but it comes at a
> > >> readability and possibly parallel cost if the global Vec needs to be
> > >> communicated to local vectors.)
> > >>
> > >>
> > > design paramaters are fixed troughout an adjoint/tlm run. They can be
> > > communicated locally once at the beginning of the run.
> > > This is what TSSetUpFromDesign and TSSetSetUpFromDesign are supposed to
> > > handle, if I get your comment.
> >
> > My point is that users currently get design parameters out of the
> > context when evaluating their RHSFunction and friends.  If that is the
> > endorsed way to access design variables, then your new function doesn't
> > need to pass the vector.  If you need to pass the parameter vector in
> > that one function, instead of obtaining them from the context, then
> > you'd need to pass the parameter vector everywhere and discourage using
> > the context for active design variables.  I think there are merits to
> > both approaches, but it absolutely needs to be consistent.
> >
> > So, to be consistent, we have to force users to perform an operation in a 
> > single way?
> > Yes, TSSetUpFromDesign is among the last things I have added, and allows to 
> > update the application context (among other things, as it is very general). 
> > I can remove the parameter vector from TSEvalGradientDAE and 
> > TSEvalGradientIC; however, having these vectors there makes it clear that 
> > we allow non-linear dependency on the parameters too. I can add a comment 
> > on the man pages that the vectors are guaranteed to be the same one passed 
> > in my TSSetUpFromDesign, or remove them. Your call.
> >
> > Users can do anything they want with the forward model context, but they 
> > are not free to change the application context of the adjoints TS. Maybe 
> > this should be improved by adding and extra slot to AdjointTSCtx to carry 
> > over the  user context (for the adjoint I mean)?
> >
> >
> > > https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d46cf7cad52/src/ts/interface/tspdeconstrainedutils.c?at=stefano_zampini%2Ffeature-continuousadjoint=file-view-default#tspdeconstrainedutils.c-579
> > >
> > > Here is how ex23.c uses it
> > >
> > > 

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Stefano Zampini]

>
>Yes, we are waiting for Hong to provide his proposed API.
>
>BTW: Stefano has not really submitted a proposed API as I requested,
> just some minor explanation of his current API and minor changes.
>
>
My proposal would be based on what we are currently discussing, and it will
build on the existing one.



-- 
Stefano

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Barry Smith]

> On Oct 17, 2017, at 8:19 AM, Jed Brown  wrote:
> 
> Matthew Knepley  writes:
> 
>> On Tue, Oct 17, 2017 at 8:51 AM, Jed Brown  wrote:
>> 
>>> Matthew Knepley  writes:
>>> 
> It's a recipe for confusion.  Either the parameters are never passed
> explicitly or they are always passed explicitly and should not be stored
> redundantly in the context, thus perhaps enabling some sort of higher
> level analysis that dynamically changes parameter values.  I would go
> with the former for now.
 
 
 I want to say again how much I dislike ad hoc memory layouts through
 contexts and the like. We have a perfectly good layout descriptor (DM)
 that should be used to describe data layout.
>>> 
>>> This is an independent change from the adjoint work and I think it's out
>>> of scope right now.  If we change it, it should go in its own PR.  I
>>> don't like having one PR with a smattering of non-essential changes to
>>> old interfaces bundled together with new conventions and new
>>> functionality.
>>> 
>> 
>> My understanding was that this discussion is not about a particular
>> PR, but about the interface we should have for sensitivity and optimal
>> control.
> 
> Sure, but changing the way parameters are passed isn't an essential part
> of the adjoint interface.  This particular debate spawned from Stefano's
> claim that passing the parameter vector to one of his functions was a
> feature of his interface relative to Hong's.  I replied that it was an
> inconsistency -- we should either change all interfaces to pass
> parameter vectors explicitly or we never pass them.
> 
> But any such change would be orthogonal to the merits of TSAdjointSolve
> versus TSCreateAdjointTS.  Indeed, discussion of how parameters are
> accessed is derailing the present discussion from this crucial
> difference in interfaces.  Let's cut it out of the current thread.
> 
> I want to hear rationale for TSAdjointSolve() or a plan for refactoring
> those algorithms into TSCreateAdjointTS().

   Yes, we are waiting for Hong to provide his proposed API.

   BTW: Stefano has not really submitted a proposed API as I requested, just 
some minor explanation of his current API and minor changes.

  Barry

> 
>>> Putting the parameters in a vector would enable finite differencing of
>>> the RHSFunction to obtain its dependence on parameters.  That might have
>>> high (non-scalable in number of parameters) cost, but would be less
>>> expensive than finite differencing an entire model.  Consider the
>>> scenario of 100 parameters and 1e6 state variables (at each time step).
>>> If we have the ability to apply the transpose of the Jacobian with
>>> respect to model state, we could run an adjoint simulation and only need
>>> 100 RHSFunction evaluations per stage, rather than 100 solves.
>>> 
>> 
>> I am not sure of the point of the above paragraph. Saying the point rather
>> than implying
>> the point helps me.
> 
> It is a potential functional case for making the parameter vector an
> explicit argument.

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Stefano Zampini]

>
> Sure, but changing the way parameters are passed isn't an essential part
> of the adjoint interface.  This particular debate spawned from Stefano's
> claim that passing the parameter vector to one of his functions was a
> feature of his interface relative to Hong's.  I replied that it was an
> inconsistency -- we should either change all interfaces to pass
> parameter vectors explicitly or we never pass them.
>
>
They could be removed from the callbacks for
TSSet{Gradient|Hessian}{DAE|IC}; however, I'd prefer to keep them in the
objective evaluation routines from clarity.




-- 
Stefano

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Matthew Knepley  writes:

> On Tue, Oct 17, 2017 at 8:51 AM, Jed Brown  wrote:
>
>> Matthew Knepley  writes:
>>
>> >> It's a recipe for confusion.  Either the parameters are never passed
>> >> explicitly or they are always passed explicitly and should not be stored
>> >> redundantly in the context, thus perhaps enabling some sort of higher
>> >> level analysis that dynamically changes parameter values.  I would go
>> >> with the former for now.
>> >
>> >
>> > I want to say again how much I dislike ad hoc memory layouts through
>> > contexts and the like. We have a perfectly good layout descriptor (DM)
>> > that should be used to describe data layout.
>>
>> This is an independent change from the adjoint work and I think it's out
>> of scope right now.  If we change it, it should go in its own PR.  I
>> don't like having one PR with a smattering of non-essential changes to
>> old interfaces bundled together with new conventions and new
>> functionality.
>>
>
> My understanding was that this discussion is not about a particular
> PR, but about the interface we should have for sensitivity and optimal
> control.

Sure, but changing the way parameters are passed isn't an essential part
of the adjoint interface.  This particular debate spawned from Stefano's
claim that passing the parameter vector to one of his functions was a
feature of his interface relative to Hong's.  I replied that it was an
inconsistency -- we should either change all interfaces to pass
parameter vectors explicitly or we never pass them.

But any such change would be orthogonal to the merits of TSAdjointSolve
versus TSCreateAdjointTS.  Indeed, discussion of how parameters are
accessed is derailing the present discussion from this crucial
difference in interfaces.  Let's cut it out of the current thread.

I want to hear rationale for TSAdjointSolve() or a plan for refactoring
those algorithms into TSCreateAdjointTS().

>> Putting the parameters in a vector would enable finite differencing of
>> the RHSFunction to obtain its dependence on parameters.  That might have
>> high (non-scalable in number of parameters) cost, but would be less
>> expensive than finite differencing an entire model.  Consider the
>> scenario of 100 parameters and 1e6 state variables (at each time step).
>> If we have the ability to apply the transpose of the Jacobian with
>> respect to model state, we could run an adjoint simulation and only need
>> 100 RHSFunction evaluations per stage, rather than 100 solves.
>>
>
> I am not sure of the point of the above paragraph. Saying the point rather
> than implying
> the point helps me.

It is a potential functional case for making the parameter vector an
explicit argument.

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Matthew Knepley]

On Tue, Oct 17, 2017 at 8:51 AM, Jed Brown  wrote:

> Matthew Knepley  writes:
>
> >> It's a recipe for confusion.  Either the parameters are never passed
> >> explicitly or they are always passed explicitly and should not be stored
> >> redundantly in the context, thus perhaps enabling some sort of higher
> >> level analysis that dynamically changes parameter values.  I would go
> >> with the former for now.
> >
> >
> > I want to say again how much I dislike ad hoc memory layouts through
> > contexts and the like. We have a perfectly good layout descriptor (DM)
> > that should be used to describe data layout.
>
> This is an independent change from the adjoint work and I think it's out
> of scope right now.  If we change it, it should go in its own PR.  I
> don't like having one PR with a smattering of non-essential changes to
> old interfaces bundled together with new conventions and new
> functionality.
>

My understanding was that this discussion is not about a particular PR, but
about
the interface we should have for sensitivity and optimal control.


> Putting the parameters in a vector would enable finite differencing of
> the RHSFunction to obtain its dependence on parameters.  That might have
> high (non-scalable in number of parameters) cost, but would be less
> expensive than finite differencing an entire model.  Consider the
> scenario of 100 parameters and 1e6 state variables (at each time step).
> If we have the ability to apply the transpose of the Jacobian with
> respect to model state, we could run an adjoint simulation and only need
> 100 RHSFunction evaluations per stage, rather than 100 solves.
>

I am not sure of the point of the above paragraph. Saying the point rather
than implying
the point helps me.

   Matt

-- 
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/ 

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Stefano Zampini]

>What happens if you call TSCreateAdjointsTS() on a TS obtained with
> TSCreateAdjointsTS()? Is the resulting TS useful?
>
>
I don't think so. Theoretically, it should be a tangent linear model with a
different forcing function




> Barry
>
>
> > On Oct 17, 2017, at 12:37 AM, Stefano Zampini 
> wrote:
> >
> >
> >
> >
> > >
> > >
> > >> In case of multiple objectives, there may be a performance reason to
> > >> amortize evaluation of several at once, though the list interface is
> > >> convenient.  Consider common objectives being quantities like lift and
> > >> drag on different surfaces of a fluids simulation or stress/strain at
> > >> certain critical joints in a structure.  Although these have some
> > >> locality, it's reasonable to assume that state dependence will have
> > >> quickly become global, thus make no attempt to handle sparse
> > >> representations of the adjoint vectors lambda.
> > >>
> > >>
> > > I don't get this comment. Is it related with multi-objective
> optimization
> > > (e.g. Pareto)?
> >
> > Adjoints are usually preferred any time you are differentiating a small
> > number of output variables with respect to a large number of inputs.  It
> > could be for multi-objective optimization, but it's every bit as
> > relevant for physical sensitivities.
> >
> >
> > If we are not talking about Pareto optimization, and thus we don't need
> a separate output from each function, then users can pass a single function
> that computes all the quantities they need at the same time. Anyway, I
> don't mind having a single callback for multiple functions.
> >
> > >> How are parameters accessed in TSComputeRHSFunction?  It looks like
> > >> they're coming out of the context.  Why should this be different?  (If
> > >> parameters need to go into a Vec, we could do that, but it comes at a
> > >> readability and possibly parallel cost if the global Vec needs to be
> > >> communicated to local vectors.)
> > >>
> > >>
> > > design paramaters are fixed troughout an adjoint/tlm run. They can be
> > > communicated locally once at the beginning of the run.
> > > This is what TSSetUpFromDesign and TSSetSetUpFromDesign are supposed to
> > > handle, if I get your comment.
> >
> > My point is that users currently get design parameters out of the
> > context when evaluating their RHSFunction and friends.  If that is the
> > endorsed way to access design variables, then your new function doesn't
> > need to pass the vector.  If you need to pass the parameter vector in
> > that one function, instead of obtaining them from the context, then
> > you'd need to pass the parameter vector everywhere and discourage using
> > the context for active design variables.  I think there are merits to
> > both approaches, but it absolutely needs to be consistent.
> >
> > So, to be consistent, we have to force users to perform an operation in
> a single way?
> > Yes, TSSetUpFromDesign is among the last things I have added, and allows
> to update the application context (among other things, as it is very
> general). I can remove the parameter vector from TSEvalGradientDAE and
> TSEvalGradientIC; however, having these vectors there makes it clear that
> we allow non-linear dependency on the parameters too. I can add a comment
> on the man pages that the vectors are guaranteed to be the same one passed
> in my TSSetUpFromDesign, or remove them. Your call.
> >
> > Users can do anything they want with the forward model context, but they
> are not free to change the application context of the adjoints TS. Maybe
> this should be improved by adding and extra slot to AdjointTSCtx to carry
> over the  user context (for the adjoint I mean)?
> >
> >
> > > https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d
> 46cf7cad52/src/ts/interface/tspdeconstrainedutils.c?at=
> stefano_zampini%2Ffeature-continuousadjoint=file-view-default#
> tspdeconstrainedutils.c-579
> > >
> > > Here is how ex23.c uses it
> > >
> > > https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d
> 46cf7cad52/src/ts/examples/tutorials/ex23.c?at=stefano_zampini%2Ffeature-
> continuousadjoint=file-view-default#ex23.c-677
> >
> > And yet you redo the scatter here instead of using what you stuffed into
> > the context.  If you needed to redo it for correctness, you'd also need
> > to in every other function that accesses design parameters.
> > https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d
> 46cf7cad52/src/ts/examples/tutorials/ex23.c?at=stefano_zampini%2Ffeature-
> continuousadjoint=file-view-default#ex23.c-274
> >
> >
> > This is a leftover from a previous version of the code (there's also a
> comment) that was not using TSSetSetUpFromDesign and it's definitely not
> needed.
> >
> >
> > >> > Both methods need the Jacobian of the DAE wrt the parameters: H
> > >> > TSAdjointSetRHSJacobian(), S TSSetGradientDAE()
> > >> >
> > >> > Initial condition dependence on the parameters is implicitly
> computed 

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Stefano Zampini]

>
>
> I don't follow why multi-objective optimization would be any different.
> Are you supposing that the different objective functions would be
> implemented by different software packages?
>
>
It was related on how I was thinking to support multi-objective
optimization (and thus multi-gradient) with the current AdjointTS, since
each adjoint run would definitely have a different adapted time step.


>
> There are a bunch of other places that allow nonlinear dependence on
> parameters yet don't pass that parameter Vec, including
> TSSetRHSFunction.  If it is passed explicitly in one such function, it
> needs to be passed explicitly in all.
>
>
Yes but TSSetRHSFunction is supposed to work for standard ODE solves (not
with the parameter dependent analysis enabled). What if instead we add a
design Vec to TS data?


> > I can add a comment on the man pages that the vectors are guaranteed
> > to be the same one passed in my TSSetUpFromDesign, or remove
> > them. Your call.
>
> It's a recipe for confusion.  Either the parameters are never passed
> explicitly or they are always passed explicitly and should not be stored
> redundantly in the context, thus perhaps enabling some sort of higher
> level analysis that dynamically changes parameter values.  I would go
> with the former for now.
>
> > Users can do anything they want with the forward model context, but they
> > are not free to change the application context of the adjoints TS. Maybe
> > this should be improved by adding and extra slot to AdjointTSCtx to carry
> > over the  user context (for the adjoint I mean)?
>
> Why do you need to own the application context for the adjoint TS,
> versus using the context parameters to TSSetRHSJacobian, etc.
>
>
You're right, I could have used the contexts. This can be fixed.


>
> Comments that some lines of code can be deleted is a good sign that they
> should be deleted.
>
>
Will do

>
> I was responding to your statement that "Initial condition dependence on
> the parameters is limited to linear dependence on all the variables" in
> his approach.  General dependence is supported, the user is just
> expected to call it themselves rather than registering a callback that
> TS calls.
>
>
This implies some knowledge from the users about how to take this gradient,
which you cannot assume in general.
If you can automatize it, why not?



>
>
> I see reasons for the adjoint TS to hold a reference to the forward TS,
> but not vice-versa.  What is the use case for your proposal below?
>
> > TSGetAdjointTS(TS f,TS* a) {
> > if (!f->adjtts) TSCreateAdjointsTS(f,>adjts);
> > *a = f->adjts:
> > }
> >
> > and the corresponding Setter to allow users to do
> >
> > TSCreateAdjointTS(ts,)
> > TSSetRHSJacobian(ats,...)
> > TSSetAdjointTS(ts,ats)
> >
> > or
> >
> > TSGetAdjointTS(ts,)
> > TSSetRHSJacobian(ats,...)
>

Just automating TSComputeObjectiveAndGradient() for PDAEs with user defined
adjoints,




-- 
Stefano

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Matthew Knepley]

On Tue, Oct 17, 2017 at 8:16 AM, Jed Brown  wrote:

> Stefano Zampini  writes:
>
> >> >
> >> >
> >> >> In case of multiple objectives, there may be a performance reason to
> >> >> amortize evaluation of several at once, though the list interface is
> >> >> convenient.  Consider common objectives being quantities like lift
> and
> >> >> drag on different surfaces of a fluids simulation or stress/strain at
> >> >> certain critical joints in a structure.  Although these have some
> >> >> locality, it's reasonable to assume that state dependence will have
> >> >> quickly become global, thus make no attempt to handle sparse
> >> >> representations of the adjoint vectors lambda.
> >> >>
> >> >>
> >> > I don't get this comment. Is it related with multi-objective
> optimization
> >> > (e.g. Pareto)?
> >>
> >> Adjoints are usually preferred any time you are differentiating a small
> >> number of output variables with respect to a large number of inputs.  It
> >> could be for multi-objective optimization, but it's every bit as
> >> relevant for physical sensitivities.
> >>
> >>
> > If we are not talking about Pareto optimization, and thus we don't need a
> > separate output from each function, then users can pass a single function
> > that computes all the quantities they need at the same time. Anyway, I
> > don't mind having a single callback for multiple functions.
>
> I don't follow why multi-objective optimization would be any different.
> Are you supposing that the different objective functions would be
> implemented by different software packages?
>
> >> >> How are parameters accessed in TSComputeRHSFunction?  It looks like
> >> >> they're coming out of the context.  Why should this be different?
> (If
> >> >> parameters need to go into a Vec, we could do that, but it comes at a
> >> >> readability and possibly parallel cost if the global Vec needs to be
> >> >> communicated to local vectors.)
> >> >>
> >> >>
> >> > design paramaters are fixed troughout an adjoint/tlm run. They can be
> >> > communicated locally once at the beginning of the run.
> >> > This is what TSSetUpFromDesign and TSSetSetUpFromDesign are supposed
> to
> >> > handle, if I get your comment.
> >>
> >> My point is that users currently get design parameters out of the
> >> context when evaluating their RHSFunction and friends.  If that is the
> >> endorsed way to access design variables, then your new function doesn't
> >> need to pass the vector.  If you need to pass the parameter vector in
> >> that one function, instead of obtaining them from the context, then
> >> you'd need to pass the parameter vector everywhere and discourage using
> >> the context for active design variables.  I think there are merits to
> >> both approaches, but it absolutely needs to be consistent.
> >>
> >
> > So, to be consistent, we have to force users to perform an operation in a
> > single way?
> > Yes, TSSetUpFromDesign is among the last things I have added, and allows
> to
> > update the application context (among other things, as it is very
> general).
> > I can remove the parameter vector from TSEvalGradientDAE and
> > TSEvalGradientIC; however, having these vectors there makes it clear that
> > we allow non-linear dependency on the parameters too.
>
> There are a bunch of other places that allow nonlinear dependence on
> parameters yet don't pass that parameter Vec, including
> TSSetRHSFunction.  If it is passed explicitly in one such function, it
> needs to be passed explicitly in all.
>
> > I can add a comment on the man pages that the vectors are guaranteed
> > to be the same one passed in my TSSetUpFromDesign, or remove
> > them. Your call.
>
> It's a recipe for confusion.  Either the parameters are never passed
> explicitly or they are always passed explicitly and should not be stored
> redundantly in the context, thus perhaps enabling some sort of higher
> level analysis that dynamically changes parameter values.  I would go
> with the former for now.


I want to say again how much I dislike ad hoc memory layouts through
contexts
and the like. We have a perfectly good layout descriptor (DM) that should
be used
to describe data layout.

  Matt


>
> > Users can do anything they want with the forward model context, but they
> > are not free to change the application context of the adjoints TS. Maybe
> > this should be improved by adding and extra slot to AdjointTSCtx to carry
> > over the  user context (for the adjoint I mean)?
>
> Why do you need to own the application context for the adjoint TS,
> versus using the context parameters to TSSetRHSJacobian, etc.
>
> >> > https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d
> >> 46cf7cad52/src/ts/interface/tspdeconstrainedutils.c?at=
> >> stefano_zampini%2Ffeature-continuousadjoint=
> file-view-default#
> >> tspdeconstrainedutils.c-579
> >> >
> >> > Here is how ex23.c uses it
> >> >
> >> > 

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Barry Smith]

  Stefano,

   What happens if you call TSCreateAdjointsTS() on a TS obtained with 
TSCreateAdjointsTS()? Is the resulting TS useful?

Barry


> On Oct 17, 2017, at 12:37 AM, Stefano Zampini  
> wrote:
> 
> 
> 
> 
> >
> >
> >> In case of multiple objectives, there may be a performance reason to
> >> amortize evaluation of several at once, though the list interface is
> >> convenient.  Consider common objectives being quantities like lift and
> >> drag on different surfaces of a fluids simulation or stress/strain at
> >> certain critical joints in a structure.  Although these have some
> >> locality, it's reasonable to assume that state dependence will have
> >> quickly become global, thus make no attempt to handle sparse
> >> representations of the adjoint vectors lambda.
> >>
> >>
> > I don't get this comment. Is it related with multi-objective optimization
> > (e.g. Pareto)?
> 
> Adjoints are usually preferred any time you are differentiating a small
> number of output variables with respect to a large number of inputs.  It
> could be for multi-objective optimization, but it's every bit as
> relevant for physical sensitivities.
> 
> 
> If we are not talking about Pareto optimization, and thus we don't need a 
> separate output from each function, then users can pass a single function 
> that computes all the quantities they need at the same time. Anyway, I don't 
> mind having a single callback for multiple functions.
>  
> >> How are parameters accessed in TSComputeRHSFunction?  It looks like
> >> they're coming out of the context.  Why should this be different?  (If
> >> parameters need to go into a Vec, we could do that, but it comes at a
> >> readability and possibly parallel cost if the global Vec needs to be
> >> communicated to local vectors.)
> >>
> >>
> > design paramaters are fixed troughout an adjoint/tlm run. They can be
> > communicated locally once at the beginning of the run.
> > This is what TSSetUpFromDesign and TSSetSetUpFromDesign are supposed to
> > handle, if I get your comment.
> 
> My point is that users currently get design parameters out of the
> context when evaluating their RHSFunction and friends.  If that is the
> endorsed way to access design variables, then your new function doesn't
> need to pass the vector.  If you need to pass the parameter vector in
> that one function, instead of obtaining them from the context, then
> you'd need to pass the parameter vector everywhere and discourage using
> the context for active design variables.  I think there are merits to
> both approaches, but it absolutely needs to be consistent.
> 
> So, to be consistent, we have to force users to perform an operation in a 
> single way? 
> Yes, TSSetUpFromDesign is among the last things I have added, and allows to 
> update the application context (among other things, as it is very general). I 
> can remove the parameter vector from TSEvalGradientDAE and TSEvalGradientIC; 
> however, having these vectors there makes it clear that we allow non-linear 
> dependency on the parameters too. I can add a comment on the man pages that 
> the vectors are guaranteed to be the same one passed in my TSSetUpFromDesign, 
> or remove them. Your call.
> 
> Users can do anything they want with the forward model context, but they are 
> not free to change the application context of the adjoints TS. Maybe this 
> should be improved by adding and extra slot to AdjointTSCtx to carry over the 
>  user context (for the adjoint I mean)?
> 
> 
> > https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d46cf7cad52/src/ts/interface/tspdeconstrainedutils.c?at=stefano_zampini%2Ffeature-continuousadjoint=file-view-default#tspdeconstrainedutils.c-579
> >
> > Here is how ex23.c uses it
> >
> > https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d46cf7cad52/src/ts/examples/tutorials/ex23.c?at=stefano_zampini%2Ffeature-continuousadjoint=file-view-default#ex23.c-677
> 
> And yet you redo the scatter here instead of using what you stuffed into
> the context.  If you needed to redo it for correctness, you'd also need
> to in every other function that accesses design parameters.
> https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d46cf7cad52/src/ts/examples/tutorials/ex23.c?at=stefano_zampini%2Ffeature-continuousadjoint=file-view-default#ex23.c-274
> 
> 
> This is a leftover from a previous version of the code (there's also a 
> comment) that was not using TSSetSetUpFromDesign and it's definitely not 
> needed. 
>  
>  
> >> > Both methods need the Jacobian of the DAE wrt the parameters: H
> >> > TSAdjointSetRHSJacobian(), S TSSetGradientDAE()
> >> >
> >> > Initial condition dependence on the parameters is implicitly computed in
> >> > Hong's code (limited to linear dependence on all the variables);
> >>
> >> How so?  Once the user gets \lambda(time=0), they can apply the chain
> >> rule to produce any dependency on the parameter vector?
> >>
> >>  Yes, the 

Re: [petsc-dev] Discussion about time-dependent optimization moved from PR

[Jed Brown]

Stefano Zampini  writes:

>> >
>> >
>> >> In case of multiple objectives, there may be a performance reason to
>> >> amortize evaluation of several at once, though the list interface is
>> >> convenient.  Consider common objectives being quantities like lift and
>> >> drag on different surfaces of a fluids simulation or stress/strain at
>> >> certain critical joints in a structure.  Although these have some
>> >> locality, it's reasonable to assume that state dependence will have
>> >> quickly become global, thus make no attempt to handle sparse
>> >> representations of the adjoint vectors lambda.
>> >>
>> >>
>> > I don't get this comment. Is it related with multi-objective optimization
>> > (e.g. Pareto)?
>>
>> Adjoints are usually preferred any time you are differentiating a small
>> number of output variables with respect to a large number of inputs.  It
>> could be for multi-objective optimization, but it's every bit as
>> relevant for physical sensitivities.
>>
>>
> If we are not talking about Pareto optimization, and thus we don't need a
> separate output from each function, then users can pass a single function
> that computes all the quantities they need at the same time. Anyway, I
> don't mind having a single callback for multiple functions.

I don't follow why multi-objective optimization would be any different.
Are you supposing that the different objective functions would be
implemented by different software packages?

>> >> How are parameters accessed in TSComputeRHSFunction?  It looks like
>> >> they're coming out of the context.  Why should this be different?  (If
>> >> parameters need to go into a Vec, we could do that, but it comes at a
>> >> readability and possibly parallel cost if the global Vec needs to be
>> >> communicated to local vectors.)
>> >>
>> >>
>> > design paramaters are fixed troughout an adjoint/tlm run. They can be
>> > communicated locally once at the beginning of the run.
>> > This is what TSSetUpFromDesign and TSSetSetUpFromDesign are supposed to
>> > handle, if I get your comment.
>>
>> My point is that users currently get design parameters out of the
>> context when evaluating their RHSFunction and friends.  If that is the
>> endorsed way to access design variables, then your new function doesn't
>> need to pass the vector.  If you need to pass the parameter vector in
>> that one function, instead of obtaining them from the context, then
>> you'd need to pass the parameter vector everywhere and discourage using
>> the context for active design variables.  I think there are merits to
>> both approaches, but it absolutely needs to be consistent.
>>
>
> So, to be consistent, we have to force users to perform an operation in a
> single way?
> Yes, TSSetUpFromDesign is among the last things I have added, and allows to
> update the application context (among other things, as it is very general).
> I can remove the parameter vector from TSEvalGradientDAE and
> TSEvalGradientIC; however, having these vectors there makes it clear that
> we allow non-linear dependency on the parameters too. 

There are a bunch of other places that allow nonlinear dependence on
parameters yet don't pass that parameter Vec, including
TSSetRHSFunction.  If it is passed explicitly in one such function, it
needs to be passed explicitly in all.

> I can add a comment on the man pages that the vectors are guaranteed
> to be the same one passed in my TSSetUpFromDesign, or remove
> them. Your call.

It's a recipe for confusion.  Either the parameters are never passed
explicitly or they are always passed explicitly and should not be stored
redundantly in the context, thus perhaps enabling some sort of higher
level analysis that dynamically changes parameter values.  I would go
with the former for now.

> Users can do anything they want with the forward model context, but they
> are not free to change the application context of the adjoints TS. Maybe
> this should be improved by adding and extra slot to AdjointTSCtx to carry
> over the  user context (for the adjoint I mean)?

Why do you need to own the application context for the adjoint TS,
versus using the context parameters to TSSetRHSJacobian, etc.

>> > https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d
>> 46cf7cad52/src/ts/interface/tspdeconstrainedutils.c?at=
>> stefano_zampini%2Ffeature-continuousadjoint=file-view-default#
>> tspdeconstrainedutils.c-579
>> >
>> > Here is how ex23.c uses it
>> >
>> > https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d
>> 46cf7cad52/src/ts/examples/tutorials/ex23.c?at=stefano_zampini%2Ffeature-
>> continuousadjoint=file-view-default#ex23.c-677
>>
>> And yet you redo the scatter here instead of using what you stuffed into
>> the context.  If you needed to redo it for correctness, you'd also need
>> to in every other function that accesses design parameters.
>>
> https://bitbucket.org/petsc/petsc/src/c2e9112e7fdfd89985f9ffc4d68b0d
>>