Re: [petsc-dev] Improving nightly builds for maint

2017-10-18 Thread Satish Balay
On Wed, 18 Oct 2017, Lisandro Dalcin wrote:

> On 18 October 2017 at 23:26, Satish Balay  wrote:
> >
> > BTW: Just so clarify - eventhough you see a 'killed' message -
> > currently those jobs are running to completion. [however long it takes
> > under valgrind]
> >
> .
> And then we get a yellow row in the nightly build matrix, right?

Pushed my change for now. [until a better fix is done]

Satish


Re: [petsc-dev] Improving nightly builds for maint

2017-10-18 Thread Lisandro Dalcin
On 18 October 2017 at 23:26, Satish Balay  wrote:
>
> BTW: Just so clarify - eventhough you see a 'killed' message -
> currently those jobs are running to completion. [however long it takes
> under valgrind]
>
.
And then we get a yellow row in the nightly build matrix, right?


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-dev] Improving nightly builds for maint

2017-10-18 Thread Satish Balay
On Wed, 18 Oct 2017, Lisandro Dalcin wrote:

> On 18 October 2017 at 22:36, Satish Balay  wrote:
> > On Wed, 18 Oct 2017, Lisandro Dalcin wrote:
> >
> > Alternate is to remove/increase the timout. It doesn't work for some
> > builds anyway [esp valgrind runs]
> >
> > This feature relies on sending 'kill' message to the process it
> > started - i.e mpiexec. But some mpiexecs when killed - don't kill the
> > child processes cleanly]
> >
> > $ git diff |cat
> > diff --git a/bin/maint/buildtest b/bin/maint/buildtest
> > index 31b5a3a028..e4bbf96db8 100755
> > --- a/bin/maint/buildtest
> > +++ b/bin/maint/buildtest
> > @@ -104,7 +104,7 @@ if ( $? ) then
> >  else
> >set VALGRIND=`echo MPIEXEC="${PETSC_DIR}/bin/petscmpiexec -valgrind"`
> >setenv PETSCVALGRIND_OPTIONS 
> > "--suppressions=${PETSC_DIR}/bin/maint/petsc-val.supp"
> > -  set TIMEOUT=720
> > +  set TIMEOUT=7200
> >  endif
> >
> 
> No, please don't!  Be green, Satish!
> 
> What's the point of running a 3D problem under valgrind for 2 hours?


BTW: Just so clarify - eventhough you see a 'killed' message -
currently those jobs are running to completion. [however long it takes
under valgrind]

Satish


> 
> I think we should somehow support
> 
> requires: !valgrind
> 
> and add special code in the test harness (maybe using some environment
> var?) to skip running these tests...
> 
> 
> 



Re: [petsc-dev] Improving nightly builds for maint

2017-10-18 Thread Lisandro Dalcin
On 18 October 2017 at 22:36, Satish Balay  wrote:
> On Wed, 18 Oct 2017, Lisandro Dalcin wrote:
>
> Alternate is to remove/increase the timout. It doesn't work for some
> builds anyway [esp valgrind runs]
>
> This feature relies on sending 'kill' message to the process it
> started - i.e mpiexec. But some mpiexecs when killed - don't kill the
> child processes cleanly]
>
> $ git diff |cat
> diff --git a/bin/maint/buildtest b/bin/maint/buildtest
> index 31b5a3a028..e4bbf96db8 100755
> --- a/bin/maint/buildtest
> +++ b/bin/maint/buildtest
> @@ -104,7 +104,7 @@ if ( $? ) then
>  else
>set VALGRIND=`echo MPIEXEC="${PETSC_DIR}/bin/petscmpiexec -valgrind"`
>setenv PETSCVALGRIND_OPTIONS 
> "--suppressions=${PETSC_DIR}/bin/maint/petsc-val.supp"
> -  set TIMEOUT=720
> +  set TIMEOUT=7200
>  endif
>

No, please don't!  Be green, Satish!

What's the point of running a 3D problem under valgrind for 2 hours?

I think we should somehow support

requires: !valgrind

and add special code in the test harness (maybe using some environment
var?) to skip running these tests...


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-dev] Improving nightly builds for maint

2017-10-18 Thread Satish Balay
On Wed, 18 Oct 2017, Lisandro Dalcin wrote:

> I opened this PR: https://bitbucket.org/petsc/petsc/pull-requests/
> 
> However, this is unlikely to fix these failures:
> 
> http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/15/examples_maint_arch-linux-pkgs-valgrind_el6.log
> 
> Satish, can we find a way to disable these tests when running under valgrind?

src/dm/examples/tutorials/ex15.c:  requires: define(PETSC_HAVE_MPIIO) 
!define(PETSC_HAVE_LIBMSMPI)

One way is to improve the test harness to support something like:

requires: !define(PETSC_ARCH="arch-linux-pkgs-valgrind")

Alternate is to remove/increase the timout. It doesn't work for some
builds anyway [esp valgrind runs]

This feature relies on sending 'kill' message to the process it
started - i.e mpiexec. But some mpiexecs when killed - don't kill the
child processes cleanly]

$ git diff |cat
diff --git a/bin/maint/buildtest b/bin/maint/buildtest
index 31b5a3a028..e4bbf96db8 100755
--- a/bin/maint/buildtest
+++ b/bin/maint/buildtest
@@ -104,7 +104,7 @@ if ( $? ) then
 else
   set VALGRIND=`echo MPIEXEC="${PETSC_DIR}/bin/petscmpiexec -valgrind"`
   setenv PETSCVALGRIND_OPTIONS 
"--suppressions=${PETSC_DIR}/bin/maint/petsc-val.supp"
-  set TIMEOUT=720
+  set TIMEOUT=7200
 endif


Satish


[petsc-dev] Improving nightly builds for maint

2017-10-18 Thread Lisandro Dalcin
I opened this PR: https://bitbucket.org/petsc/petsc/pull-requests/

However, this is unlikely to fix these failures:

http://ftp.mcs.anl.gov/pub/petsc/nightlylogs/archive/2017/10/15/examples_maint_arch-linux-pkgs-valgrind_el6.log

Satish, can we find a way to disable these tests when running under valgrind?



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


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

2017-10-18 Thread Stefano Zampini

> On Oct 18, 2017, at 7:31 PM, Barry Smith  wrote:
> 
> 
>> On Oct 18, 2017, at 4:47 AM, Stefano Zampini  
>> wrote:
>> 
>> I don't have a real API, just a sketch of a possible general framework (with 
>> holes)
>> 
>> First, I prefer having TSCreateAdjointTS(ts,) (and possibly 
>> SNESGetAdjointSNES for optimization with nonlinear time-independent 
>> problems) that allows users to replace their Jacobian evaluation routines.
> 
>   So you don't like my idea of pulling the sensitivity API out of TS and 
> putting it into a higher level object (that serves up the needed integrators 
> via PetscSensiGetTS() and PetscSensiGetAdjointTS()?  Do you think that my 
> model adds unneeded complexity to users, or what?

I agree on pulling it out, but instead of having PetscSensi, I proposed to have 
a DM that has the API with callbacks for parameter depedent evaluation of 
residuals and Jacobians (and Hessian terms).
It can also be TSGetDMDesign() and SNESGetDMDesign() instead of DMGetDMDesign 
which does not sound very well.

TSCreateAdjointTS() just gives you back a TS that is able to perform the 
adjoint run (either continuous or discrete) and shares the DMDesign, that can 
activate first-order or second order adjoints.
Both approaches have the same forcing term API (that only differ by the 
sampling time) which samples the state derivative of the objective functions 
(for gradient computations) and on some hessian terms for second-order adjoints 
(don’t know the details for discrete adjoints) The quadrature function they 
need to perform is the same, except for the input. I think DMDesign can also 
provide support for quadrature callbacks.



> 
>  Barry
> 
> 
> 
>> 
>> Then, there's the need of carrying over the design vector to TS (and SNES), 
>> and this can be accomplished by a DM
>> 
>> // User code
>> DMGetDMDesign(dm,)
>> DMSetApplicationContext(ddm,userctx);   // attach user 
>> context
>> DMDesignSetSetUp(ddm,(PetscErrorCode*)(DM,Vec)) //setup user context
>> DMDesignSetGlobalToLocal //maybe other callbacks ???
>> DMDesignSetDesignVec(ddm,currentdesign)//current design 
>> vector, can be updated by higher level routines
>> 
>> Then, within PETSc, we can do
>> 
>> TSGetDM(ts,)
>> DMGetDMTS(dm,)
>> DMGetDMDesign(dm,)
>> DMTSSetDMDesign(tdm,ddm) // without taking reference on ddm, just passing 
>> callbacks from DMDesign to private DMTS ops 
>> 
>> and for gradient based optimization and sensitivity for DAEs, we need to 
>> implement the equivalents of TSSet{Gradient|Hessian}{DAE|IC} in DMTS. 
>> DMDesign should provide an API with callbacks for computing gradients (and 
>> possibly hessian terms) of residual evaluations wrt the parameters.
>> 
>> Similarly for a parameter dependent SNES
>> 
>> SNESGetDM(snes,)
>> DMGetDMSNES(dm,)
>> DMGetDMDesign(dm,)
>> DMSNESSetDMDesign(kdm,ddm) // without taking reference on ddm, just passing 
>> callbacks from DMDesign to private DMSNES ops
>> 
>> Then, we need a consistent way of expressing objective functions 
>> (DMTAOObjective???), and embed quadratures inside TSSolve(). Thoughts?
>> 
>> 
>> 
>> 2017-10-17 22:02 GMT+03:00 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 

Re: [petsc-dev] PetscSF in Fortran

2017-10-18 Thread Matthew Knepley
On Tue, Oct 17, 2017 at 11:35 PM, Adrian Croucher  wrote:

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


I do not understand Fortran. However, it looks like you can choose an
arbitrary type for the array, such as character,
and then use the 'transfer' intrinsic from F95 in your code to get an array
like that

  https://jblevins.org/log/transfer

so you only need an interface with that generic array type in PETSc.

   Matt


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


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

2017-10-18 Thread Barry Smith

> On Oct 18, 2017, at 4:47 AM, Stefano Zampini  
> wrote:
> 
> I don't have a real API, just a sketch of a possible general framework (with 
> holes)
> 
> First, I prefer having TSCreateAdjointTS(ts,) (and possibly 
> SNESGetAdjointSNES for optimization with nonlinear time-independent problems) 
> that allows users to replace their Jacobian evaluation routines.

   So you don't like my idea of pulling the sensitivity API out of TS and 
putting it into a higher level object (that serves up the needed integrators 
via PetscSensiGetTS() and PetscSensiGetAdjointTS()?  Do you think that my model 
adds unneeded complexity to users, or what?

  Barry



> 
> Then, there's the need of carrying over the design vector to TS (and SNES), 
> and this can be accomplished by a DM
> 
> // User code
> DMGetDMDesign(dm,)
> DMSetApplicationContext(ddm,userctx);   // attach user 
> context
> DMDesignSetSetUp(ddm,(PetscErrorCode*)(DM,Vec)) //setup user context
> DMDesignSetGlobalToLocal //maybe other callbacks ???
> DMDesignSetDesignVec(ddm,currentdesign)//current design 
> vector, can be updated by higher level routines
> 
> Then, within PETSc, we can do
> 
> TSGetDM(ts,)
> DMGetDMTS(dm,)
> DMGetDMDesign(dm,)
> DMTSSetDMDesign(tdm,ddm) // without taking reference on ddm, just passing 
> callbacks from DMDesign to private DMTS ops 
> 
> and for gradient based optimization and sensitivity for DAEs, we need to 
> implement the equivalents of TSSet{Gradient|Hessian}{DAE|IC} in DMTS. 
> DMDesign should provide an API with callbacks for computing gradients (and 
> possibly hessian terms) of residual evaluations wrt the parameters.
> 
> Similarly for a parameter dependent SNES
> 
> SNESGetDM(snes,)
> DMGetDMSNES(dm,)
> DMGetDMDesign(dm,)
> DMSNESSetDMDesign(kdm,ddm) // without taking reference on ddm, just passing 
> callbacks from DMDesign to private DMSNES ops
> 
> Then, we need a consistent way of expressing objective functions 
> (DMTAOObjective???), and embed quadratures inside TSSolve(). Thoughts?
> 
> 
> 
> 2017-10-17 22:02 GMT+03:00 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 

[petsc-dev] Cleanup of dmhooks in snesdestroy (was: Re: Discussion about time-dependent optimization moved from PR)

2017-10-18 Thread Lawrence Mitchell

> On 17 Oct 2017, at 17:40, Lawrence Mitchell  wrote:
> 
> 
>> 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.

Here we go.

This just reuses a DA in a loop making two SNESes to solve the same problem.

$ ./dmhookcleanup -da_refine 1 -pc_type mg  -snes_fd_color -snes_monitor
  0 SNES Function norm 4.3750e-01 
  1 SNES Function norm 1.086628858540e-05 
  2 SNES Function norm 7.356279464858e-12 
  0 SNES Function norm 4.3750e-01 
[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Wrong type of object: Parameter # 2
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for 
trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.7.6-3851-g17f6384  GIT Date: 
2017-05-09 22:45:05 -0500
[0]PETSC ERROR: ./dmhookcleanup on a arch-linux2-c-dbg named yam.doc.ic.ac.uk 
by lmitche1 Wed Oct 18 11:04:11 2017
[0]PETSC ERROR: Configure options --download-chaco=1 --download-ctetgen=1 
--download-exodusii=1 --download-hdf5=1 --download-hypre=1 --download-metis=1 
--download-ml=1 --download-mumps=1 --download-netcdf=1 --download-parmetis=1 
--download-ptscotch=1 --download-scalapack=1 --download-triangle=1 
--download-eigen=1 --with-c2html=0 --with-debugging=1 --with-make-np=32 
--with-shared-libraries=1 PETSC_ARCH=arch-linux2-c-dbg
[0]PETSC ERROR: #1 MatRestrict() line 8044 in 
/data/lmitche1/src/deps/petsc/src/mat/interface/matrix.c
[0]PETSC ERROR: #2 DMRestrictHook_SNESVecSol() line 528 in 
/data/lmitche1/src/deps/petsc/src/snes/interface/snes.c
[0]PETSC ERROR: #3 DMRestrict() line 2421 in 
/data/lmitche1/src/deps/petsc/src/dm/interface/dm.c
[0]PETSC ERROR: #4 PCSetUp_MG() line 729 in 
/data/lmitche1/src/deps/petsc/src/ksp/pc/impls/mg/mg.c
[0]PETSC ERROR: #5 PCSetUp() line 924 in 
/data/lmitche1/src/deps/petsc/src/ksp/pc/interface/precon.c
[0]PETSC ERROR: #6 KSPSetUp() line 378 in 
/data/lmitche1/src/deps/petsc/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #7 KSPSolve() line 609 in 
/data/lmitche1/src/deps/petsc/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #8 SNESSolve_NEWTONLS() line 224 in 
/data/lmitche1/src/deps/petsc/src/snes/impls/ls/ls.c
[0]PETSC ERROR: #9 SNESSolve() line 4106 in 
/data/lmitche1/src/deps/petsc/src/snes/interface/snes.c
[0]PETSC ERROR: #10 main() line 61 in 
/data/lmitche1/src/petsc-doodles/dmhookcleanup.c
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -da_refine 1
[0]PETSC ERROR: -pc_type mg
[0]PETSC ERROR: -snes_fd_color
[0]PETSC ERROR: -snes_monitor
[0]PETSC ERROR: End of Error Message ---send entire error 
message to petsc-ma...@mcs.anl.gov--
--
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD 
with errorcode 62.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--

valgrind log:

==2215210== Memcheck, a memory error detector
==2215210== Copyright (C) 2002-2015, and GNU GPL'd, by Julian Seward et al.
==2215210== Using Valgrind-3.11.0 and LibVEX; rerun with -h for copyright info
==2215210== Command: ./dmhookcleanup -da_refine 1 -pc_type mg -snes_fd_color 
-snes_monitor
==2215210== Parent PID: 1523638
==2215210== 
==2215210== Warning: noted but unhandled ioctl 0x3001 with no 
size/direction hints.
==2215210==This could cause spurious value errors to appear.
==2215210==See README_MISSING_SYSCALL_OR_IOCTL for guidance on writing a 
proper wrapper.
==2215210== Warning: noted but unhandled ioctl 0x27 with no size/direction 
hints.
==2215210==This could cause spurious value errors to appear.
==2215210==See README_MISSING_SYSCALL_OR_IOCTL for guidance on writing a 
proper wrapper.
==2215210== Warning: noted but unhandled ioctl 0x7ff with no size/direction 
hints.
==2215210==This could cause spurious value errors to appear.
==2215210==See README_MISSING_SYSCALL_OR_IOCTL for guidance on writing a 
proper wrapper.
==2215210== Warning: noted but unhandled ioctl 0x25 with no size/direction 
hints.
==2215210==This could cause spurious value errors to appear.
==2215210==See README_MISSING_SYSCALL_OR_IOCTL for guidance on writing a 
proper wrapper.
==2215210== Warning: noted but unhandled ioctl 0x17 with no size/direction 
hints.
==2215210==This could cause spurious value errors to appear.
==2215210==See README_MISSING_SYSCALL_OR_IOCTL for 

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

2017-10-18 Thread Stefano Zampini
I don't have a real API, just a sketch of a possible general framework
(with holes)

First, I prefer having TSCreateAdjointTS(ts,) (and possibly
SNESGetAdjointSNES for optimization with nonlinear time-independent
problems) that allows users to replace their Jacobian evaluation routines.

Then, there's the need of carrying over the design vector to TS (and SNES),
and this can be accomplished by a DM

// User code
DMGetDMDesign(dm,)
DMSetApplicationContext(ddm,userctx);   // attach user
context
DMDesignSetSetUp(ddm,(PetscErrorCode*)(DM,Vec)) //setup user context
DMDesignSetGlobalToLocal //maybe other callbacks ???
DMDesignSetDesignVec(ddm,currentdesign)//current design
vector, can be updated by higher level routines

Then, within PETSc, we can do

TSGetDM(ts,)
DMGetDMTS(dm,)
DMGetDMDesign(dm,)
DMTSSetDMDesign(tdm,ddm) // without taking reference on ddm, just passing
callbacks from DMDesign to private DMTS ops

and for gradient based optimization and sensitivity for DAEs, we need to
implement the equivalents of TSSet{Gradient|Hessian}{DAE|IC} in DMTS.
DMDesign should provide an API with callbacks for computing gradients (and
possibly hessian terms) of residual evaluations wrt the parameters.

Similarly for a parameter dependent SNES

SNESGetDM(snes,)
DMGetDMSNES(dm,)
DMGetDMDesign(dm,)
DMSNESSetDMDesign(kdm,ddm) // without taking reference on ddm, just passing
callbacks from DMDesign to private DMSNES ops

Then, we need a consistent way of expressing objective functions
(DMTAOObjective???), and embed quadratures inside TSSolve(). Thoughts?



2017-10-17 22:02 GMT+03:00 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