date:20171020

Barry Smith  writes:

>> What does that has to do with PetscSF, which never used PetscDataType
>> and uses MPI datatypes with support for derived types (albeit with a
>> limited number of combiners)?
>
>Because the code he started to write uses F90Array1dAccess() which
>is built around PetscDataType, hence he needs to convert the
>MPI_Datatype. Which is where all the trouble comes from. There are
>multiple ways of resolving his problem, it more depends on my time
>then anything else since it is unlikely anyone but you or I will
>"fix" this, and you don't have time.

Ah, okay.  Could we write a type generic interface for that now that we
use more recent Fortran?  We could pass an MPI_Datatype and just check
sizes instead of trying to fully resolve the type.  I think this is
possible with F2008+TS29113, but I haven't looked at those bindings for
a few years and don't remember details.

Re: [petsc-dev] SuperLU failure with valgrind


  Mark,

   Run valgrind on a vanilla Linux system, valgrind is not great on Macs. It 
should be very useful there. And give you a line number.

I doubt it is a false positive. 

   Barry

> On Oct 16, 2017, at 8:11 PM, Mark Adams  wrote:
> 
> Yep, there was a long thread on segvs with pdgssvx in SuperLU_dist almost 
> exactly a year ago. Good catch Matt: "SuperLU_dist issue in 3.7.4".
> 
> Its not clear that it was resolved, but the code works, just barfs in 
> valgrind on my osx. And Valgrind is barfing on options data base methods on 
> Cori. 
> 
> Again the code runs fine though. Probably false positives.
> 
> On Mon, Oct 16, 2017 at 12:31 PM, Matthew Knepley  wrote:
> We had a previous error with pdgssvx in SuperLU I think. Maybe searching 
> petsc-maint would get it?
> 
>Matt
> 
> On Mon, Oct 16, 2017 at 12:21 PM, Mark Adams  wrote:
> I just ran this and have a little bit of a stack trace. This is on my laptop 
> and MPI can be a little flaky here (eg, IBarrier does not work). I am going 
> to move to Cori soon and I will try to reproduce this.
> Thanks,
> 
> ==68941==at 0x103A66AA8: MPIR_Process_status (mpiimpl.h:4394)
> ==68941==by 0x103A6852F: MPIC_Waitall (helper_fns.c:774)
> ==68941==by 0x1038ECE88: MPIR_Alltoallv_intra (alltoallv.c:194)
> ==68941==by 0x1038ED7F9: MPIR_Alltoallv (alltoallv.c:339)
> ==68941==by 0x1038EDA53: MPIR_Alltoallv_impl (alltoallv.c:376)
> ==68940==at 0x103A66AA8: MPIR_Process_status (mpiimpl.h:4394)
> ==68940==by 0x103A6852F: MPIC_Waitall (helper_fns.c:774)
> ==68940==by 0x1038ECE88: MPIR_Alltoallv_intra (alltoallv.c:194)
> ==68940==by 0x1038ED7F9: MPIR_Alltoallv (alltoallv.c:339)
> ==68940==by 0x1038EDA53: MPIR_Alltoallv_impl (alltoallv.c:376)
> ==68940==by 0x103719112: MPI_Alltoallv (alltoallv.c:527)
> ==68940==by 0x10238B87D: pdCompRow_loc_to_CompCol_global (in 
> /Users/markadams/Codes/petsc/arch-macosx-gnu-g/lib/libsuperlu_dist.5.1.3.dylib)
> ==68940==by 0x1023800CB: pdgssvx (in 
> /Users/markadams/Codes/petsc/arch-macosx-gnu-g/lib/libsuperlu_dist.5.1.3.dylib)
> ==68940==by 0x100AB92DB: MatLUFactorNumeric_SuperLU_DIST 
> (superlu_dist.c:429)
> 
> 
> On Mon, Oct 16, 2017 at 12:05 PM, Xiaoye S. Li  wrote:
> Mark,
> Is it possible to get the line number? 
> For example, the first failure is
> 
> ==63582== Conditional jump or move depends on uninitialised value(s)
> ==63582==at 0x103A5FAA8: MPIR_Process_status (mpiimpl.h:4394)
> ==63582==by 0x103A6152F: MPIC_Waitall (helper_fns.c:774)
> ==63582==by 0x1038E2A34: MPIR_Alltoall_intra (alltoall.c:369)
> ==63582==by 0x1038E35E1: MPIR_Alltoall (alltoall.c:564)
> ==63582==by 0x1038E37E6: MPIR_Alltoall_impl (alltoall.c:599)
> ==63582==by 0x1037106AD: MPI_Alltoall (alltoall.c:722)
> ==63582==by 0x10236EA7C: static_schedule (in 
> /Users/markadams/Codes/petsc/arch-macosx-gnu-g/lib/libsuperlu_dist.5.1.3.dylib)
> 
> I checked all the MPI_alltoall in static_schedule() routine, I don't see any 
> problem.
> 
> Sherry
> 
> 
> 
> On Mon, Oct 16, 2017 at 7:21 AM, Mark Adams  wrote:
> FYI, I get this error on one processor with SuperLU under valgrind. Could 
> this just be a valgrind issue?
> 
> Mark
> 
> /Users/markadams/Codes/petsc/arch-macosx-gnu-g/bin/mpiexec -n 1 valgrind 
> --dsymutil=yes --leak-check=no --gen-suppressions=no --num-callers=20 
> --error-limit=no ./ex48 -debug 2 -dim 2 -dm_refine 3 -ts_monitor -implicit 
> true -ts_type beuler -pc_type lu -pc_factor_mat_solver_package superlu_dist 
> -ksp_type preonly -snes_monitor -snes_rtol 1.e-10 -snes_stol 1.e-10 
> -snes_converged_reason -snes_atol 1.e-18 -snes_converged_reason 
> -petscspace_order 2 -petscspace_poly_tensor -ts_max_steps 1 -ts_dt 1.e-3 -eps 
> 1.e-12 -eta 0.001 -ves 0.005 -beta 0.01 -mu 0.0002 -dm_view hdf5:sol.h5 
> -vec_view hdf5:sol.h5::append -dm_plex_periodic_cut -y_periodicity PERIODIC 
> -cells 2,4 -Jop 4.99 -line_dir 1,1 -line_coord 3.14159265359,1.57079632679 
> -real_view :u.m:ascii_matlab -fft_view :spectra.m:ascii_matlab
> ==63582== Memcheck, a memory error detector
> ==63582== Copyright (C) 2002-2017, and GNU GPL'd, by Julian Seward et al.
> ==63582== Using Valgrind-3.13.0 and LibVEX; rerun with -h for copyright info
> ==63582== Command: ./ex48 -debug 2 -dim 2 -dm_refine 3 -ts_monitor -implicit 
> true -ts_type beuler -pc_type lu -pc_factor_mat_solver_package superlu_dist 
> -ksp_type preonly -snes_monitor -snes_rtol 1.e-10 -snes_stol 1.e-10 
> -snes_converged_reason -snes_atol 1.e-18 -snes_converged_reason 
> -petscspace_order 2 -petscspace_poly_tensor -ts_max_steps 1 -ts_dt 1.e-3 -eps 
> 1.e-12 -eta 0.001 -ves 0.005 -beta 0.01 -mu 0.0002 -dm_view hdf5:sol.h5 
> -vec_view hdf5:sol.h5::append -dm_plex_periodic_cut -y_periodicity PERIODIC 
> -cells 2,4 -Jop 4.99 -line_dir 1,1 -line_coord 3.14159265359,1.57079632679 
> -real_view :u.m:ascii_matlab -fft_view

Re: [petsc-dev] PetscSF in Fortran


> On Oct 20, 2017, at 8:09 PM, Jed Brown  wrote:
> 
> Barry Smith  writes:
> 
>>> On Oct 20, 2017, at 12:31 PM, Jed Brown  wrote:
>>> 
>>> Barry Smith  writes:
>>> 
  Adrian,
 
   You should not use F90Array1d *rptr as arguments in the Fortran 
 interface. You should just use a regular Fortran one dimensional array of 
 real/scalar.
 Fortran doesn't handle polymorphism in this way at all. You have to have 
 multiple f90 interfaces, one for each type and provide a C stub for real, 
 int, or whatever else you want to send.
>>> 
>>> Barry, look at the "use mpi_f08" way of calling MPI from Fortran.
>> 
>>  Jed,
>> 
>>   Rather terse response. 
>> 
>>   Are you suggesting in PETSc we use type(*) to manage multiple types 
>> through the same function?  Looks doable, I wasn't aware of this. This could 
>> possibly reduce a lot of code duplication we currently have. 
> 
> PetscSF uses MPI type handling and thus it would make sense to use a
> similarly designed Fortran module.
> 
>>   Still I would like to get rid of the use PetscDataType rather than write 
>> new code that uses it.
>> 
>>   I need to think more in this case.
>> 
>>   Waiting to hear from Adrian what types he needs to pass around (use of 
>> PetscDataType restricts to built in MPI datatypes regardless of what Fortran 
>> interface approach we use
> 
> What does that has to do with PetscSF, which never used PetscDataType
> and uses MPI datatypes with support for derived types (albeit with a
> limited number of combiners)?

   Because the code he started to write uses F90Array1dAccess() which is built 
around PetscDataType, hence he needs to convert the MPI_Datatype. Which is 
where all the trouble comes from. There are multiple ways of resolving his 
problem, it more depends on my time then anything else since it is unlikely 
anyone but you or I will "fix" this, and you don't have time.

   

  Barry

Re: [petsc-dev] PetscSF in Fortran

Barry Smith  writes:

>> On Oct 20, 2017, at 12:31 PM, Jed Brown  wrote:
>> 
>> Barry Smith  writes:
>> 
>>>   Adrian,
>>> 
>>>You should not use F90Array1d *rptr as arguments in the Fortran 
>>> interface. You should just use a regular Fortran one dimensional array of 
>>> real/scalar.
>>> Fortran doesn't handle polymorphism in this way at all. You have to have 
>>> multiple f90 interfaces, one for each type and provide a C stub for real, 
>>> int, or whatever else you want to send.
>> 
>> Barry, look at the "use mpi_f08" way of calling MPI from Fortran.
>
>   Jed,
>
>Rather terse response. 
>
>Are you suggesting in PETSc we use type(*) to manage multiple types 
> through the same function?  Looks doable, I wasn't aware of this. This could 
> possibly reduce a lot of code duplication we currently have. 

PetscSF uses MPI type handling and thus it would make sense to use a
similarly designed Fortran module.

>Still I would like to get rid of the use PetscDataType rather than write 
> new code that uses it.
>
>I need to think more in this case.
>
>Waiting to hear from Adrian what types he needs to pass around (use of 
> PetscDataType restricts to built in MPI datatypes regardless of what Fortran 
> interface approach we use

What does that has to do with PetscSF, which never used PetscDataType
and uses MPI datatypes with support for derived types (albeit with a
limited number of combiners)?

Re: [petsc-dev] PetscSF in Fortran


> On Oct 20, 2017, at 12:31 PM, Jed Brown  wrote:
> 
> Barry Smith  writes:
> 
>>   Adrian,
>> 
>>You should not use F90Array1d *rptr as arguments in the Fortran 
>> interface. You should just use a regular Fortran one dimensional array of 
>> real/scalar.
>> Fortran doesn't handle polymorphism in this way at all. You have to have 
>> multiple f90 interfaces, one for each type and provide a C stub for real, 
>> int, or whatever else you want to send.
> 
> Barry, look at the "use mpi_f08" way of calling MPI from Fortran.

  Jed,

   Rather terse response. 

   Are you suggesting in PETSc we use type(*) to manage multiple types through 
the same function?  Looks doable, I wasn't aware of this. This could possibly 
reduce a lot of code duplication we currently have. 

   Still I would like to get rid of the use PetscDataType rather than write new 
code that uses it.

   I need to think more in this case.

   Waiting to hear from Adrian what types he needs to pass around (use of 
PetscDataType restricts to built in MPI datatypes regardless of what Fortran 
interface approach we use


  Barry

> 
>>   What types do you need to send? It is probably easier if we just write the 
>> Fortran interfaces for you.
>> 
>>   If you can send a simple Fortran PETS  code that uses PetscSFBcastBegin() 
>> and PetscSFBcastEnd() that I can use for testing then I will write them. Do 
>> a bcast with int and another with real/scalar.
>> 
>> 
>>   Barry
>> 
>> 
>> 
>>> On Oct 19, 2017, at 8:59 PM, Adrian Croucher  
>>> wrote:
>>> 
>>> hi
>>> 
>>> On 19/10/17 06:45, Matthew Knepley wrote:
 On Tue, Oct 17, 2017 at 11:35 PM, Adrian Croucher 
  wrote:
 
 
 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 think this is working now by just declaring those void* C variables as 
>>> type(*) in the Fortran interface, e.g.:
>>> 
>>>  Interface
>>> Subroutine PetscSFBcastBegin(sf,unit,rarray,
>>> &   larray,ierr)
>>>   use petscisdef
>>>   PetscSF :: sf
>>>   PetscInt :: unit
>>>   type(*) :: rarray(:)
>>>   type(*) :: larray(:)
>>>   PetscErrorCode :: ierr
>>> End Subroutine PetscSFBcastBegin
>>>  End Interface
>>> 
>>> The only difficulty I have left with this is in the MPI_Datatype variable. 
>>> I'd forgotten that these datatypes are all different in C and Fortran as 
>>> well.
>>> 
>>> I amended the C interface code to the following, to convert the Fortran MPI 
>>> datatype (actually an integer) to a C MPI_Datatype:
>>> 
>>> PETSC_EXTERN void PETSC_STDCALL petscsfbcastbegin_(PetscSF *sf, MPI_Fint 
>>> *unit, F90Array1d *rptr, F90Array1d *lptr , int *ierr 
>>> PETSC_F90_2PTR_PROTO(rptrd) PETSC_F90_2PTR_PROTO(lptrd))
>>> {
>>>  PetscDataType pdtype;
>>>  MPI_Datatype dtype;
>>>  const void* rootdata;
>>>  void* leafdata;
>>> 
>>>  dtype = MPI_Type_f2c(*unit);
>>>  *ierr = PetscMPIDataTypeToPetscDataType(dtype, );if (*ierr) return;
>>>  *ierr = F90Array1dAccess(rptr, pdtype, (void**)  
>>> PETSC_F90_2PTR_PARAM(rptrd));if (*ierr) return;
>>>  *ierr = F90Array1dAccess(lptr, pdtype, (void**)  
>>> PETSC_F90_2PTR_PARAM(lptrd));if (*ierr) return;
>>> 
>>>  *ierr = PetscSFBcastBegin(*sf, dtype, rootdata, leafdata);
>>> 
>>> }
>>> 
>>> The problem is this only seems to work if I declare the datatype in the 
>>> calling Fortran code to be of the appropriate C MPI datatype, e.g. MPI_INT, 
>>> rather than the corresponding Fortran datatype, e.g. MPI_INTEGER (which 
>>> causes PetscMPIDataTypeToPetscDataType() to fail, as something weird

Re: [petsc-dev] TS Terminology

On Fri, Oct 20, 2017 at 6:44 PM, Jed Brown  wrote:

> Matthew Knepley  writes:
>
> > On Fri, Oct 20, 2017 at 6:36 PM, Jed Brown  wrote:
> >
> >> Matthew Knepley  writes:
> >>
> >> > I think you are right that the SNES thing needs fixed as well. The
> >> > problem has not come up because no one ever uses that interface. It is
> >> > only there because it makes FAS a lot easier.
> >>
> >> That interface precedes FAS.  It's there because it's convenient.
> >>
> >
> > That is not my memory.
>
> You made the change.  That was long before FAS moved from DMMG into
> SNES.
>

I know. I was writing FAS, code that no longer exists, and needed that to
make things nice.

   Matt


> commit f69a0ea3504314d029ee423e0de2950ece2dab86
> Author: Matthew Knepley 
> Date:   Tue Mar 15 18:34:40 2005 -0600
>
> bk-changeset-1.2975.1.5
> knepley@khan.(none)|ChangeSet|20050316003440|07799
> ChangeSet
>   1.2975.1.5 05/03/15 18:34:40 knepley@khan.(none) +235 -0
>   Fixed many const arguments
>   Changed MatCreate() to only take a communicator
>- Added MatSetSizes()
>   Changed SNESSolve() to accept a constant argument
>- This should allow SNES to look almost like KSP
>   Added Options methods, SetUp(), and View() and objects which lacked
> them
>
>


-- 
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] MatAssembly and VecAssembly

Barry Smith  writes:

>> On Oct 20, 2017, at 5:33 PM, Kong, Fande  wrote:
>> 
>> 
>> 
>> On Fri, Oct 20, 2017 at 4:24 PM, Barry Smith  wrote:
>> 
>> > On Oct 20, 2017, at 12:17 PM, Kong, Fande  wrote:
>> >
>> > How about to do  a global check (MPI_Allreduce)?
>> 
>>Kill you for large number of processes.
>> 
>> > If we do not set values at all and the matrix is already assembled, we 
>> > just return without doing anything?
>> >
>> > How expensive, in the current implementation,  to call 
>> > MatAssemblyBegin/End if there are no any stashed data? Is it so cheap that 
>> > we can just ignore it?
>> 
>>It requires at least one global reduction.
>> 
>> So, this global reduction is cheaper than MPI_Allreduce?  Or they are 
>> similar.
>
>   No it is not cheaper, it is the same thing. 

The point-to-point communication when using -vec_assembly_bts with
VEC_SUBSET_OFF_PROC_ENTRIES may be slightly less expensive.

>   You need to eliminate unneeded AssemblyBegin/End in MOOSE. That is the only 
> corrective action that can take place.

You should definitely do this rather than trying to hide behind possibly
slightly faster null assembly.

Re: [petsc-dev] TS Terminology

Matthew Knepley  writes:

> On Fri, Oct 20, 2017 at 6:36 PM, Jed Brown  wrote:
>
>> Matthew Knepley  writes:
>>
>> > I think you are right that the SNES thing needs fixed as well. The
>> > problem has not come up because no one ever uses that interface. It is
>> > only there because it makes FAS a lot easier.
>>
>> That interface precedes FAS.  It's there because it's convenient.
>>
>
> That is not my memory.

You made the change.  That was long before FAS moved from DMMG into
SNES.

commit f69a0ea3504314d029ee423e0de2950ece2dab86
Author: Matthew Knepley 
Date:   Tue Mar 15 18:34:40 2005 -0600

bk-changeset-1.2975.1.5
knepley@khan.(none)|ChangeSet|20050316003440|07799
ChangeSet
  1.2975.1.5 05/03/15 18:34:40 knepley@khan.(none) +235 -0
  Fixed many const arguments
  Changed MatCreate() to only take a communicator
   - Added MatSetSizes()
  Changed SNESSolve() to accept a constant argument
   - This should allow SNES to look almost like KSP
  Added Options methods, SetUp(), and View() and objects which lacked them

Re: [petsc-dev] TS Terminology


> On Oct 20, 2017, at 5:36 PM, Jed Brown  wrote:
> 
> Matthew Knepley  writes:
> 
>> I think you are right that the SNES thing needs fixed as well. The
>> problem has not come up because no one ever uses that interface. It is
>> only there because it makes FAS a lot easier.
> 
> That interface precedes FAS.  It's there because it's convenient.

  It doesn't matter when it started.

Re: [petsc-dev] TS Terminology

On Fri, Oct 20, 2017 at 6:36 PM, Jed Brown  wrote:

> Matthew Knepley  writes:
>
> > I think you are right that the SNES thing needs fixed as well. The
> > problem has not come up because no one ever uses that interface. It is
> > only there because it makes FAS a lot easier.
>
> That interface precedes FAS.  It's there because it's convenient.
>

That is not my memory.

  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] TS Terminology

Matthew Knepley  writes:

> I think you are right that the SNES thing needs fixed as well. The
> problem has not come up because no one ever uses that interface. It is
> only there because it makes FAS a lot easier.

That interface precedes FAS.  It's there because it's convenient.

Re: [petsc-dev] MatAssembly and VecAssembly


> On Oct 20, 2017, at 5:33 PM, Kong, Fande  wrote:
> 
> 
> 
> On Fri, Oct 20, 2017 at 4:24 PM, Barry Smith  wrote:
> 
> > On Oct 20, 2017, at 12:17 PM, Kong, Fande  wrote:
> >
> > How about to do  a global check (MPI_Allreduce)?
> 
>Kill you for large number of processes.
> 
> > If we do not set values at all and the matrix is already assembled, we just 
> > return without doing anything?
> >
> > How expensive, in the current implementation,  to call MatAssemblyBegin/End 
> > if there are no any stashed data? Is it so cheap that we can just ignore it?
> 
>It requires at least one global reduction.
> 
> So, this global reduction is cheaper than MPI_Allreduce?  Or they are similar.

  No it is not cheaper, it is the same thing. 

  You need to eliminate unneeded AssemblyBegin/End in MOOSE. That is the only 
corrective action that can take place.


> 
> 
>  
> >
> > I am asking because we call MatAssemblyBegin/End so often in MOOSE. I want 
> > to make sure this is not going to bring up any performance issue.
> 
>   You need to check each use in MOOSE and see WHY it is being called. If no 
> reason then don't call.
> 
> Good point.
> 
> Fande,
> 
>  
> 
>   Barry
> 
> >
> >
> > Fande,
> >
> > On Fri, Oct 20, 2017 at 11:08 AM, Barry Smith  wrote:
> >
> >   One process sets a value in the matrix, the others do not. They all call 
> > MatAssemblyBegin(). Some processes will skip the assembly and hence the 
> > code will hang.
> >
> >
> > > On Oct 20, 2017, at 12:03 PM, Kong, Fande  wrote:
> > >
> > > Hi All,
> > >
> > > In Mat/VecAssemblyBegin/End, why we do not check if or not mat/vec is 
> > > assembled. If mat/vec is already assembled, should we just return without 
> > > doing anything?
> > >
> > > I think we have some particular reasons not to check if the matrix is 
> > > assembled. I honestly do not know why.
> > >
> > > Thanks,
> > >
> > > Fande,
> >
> >

Re: [petsc-dev] MatAssembly and VecAssembly

2017-10-20 Thread Kong, Fande

On Fri, Oct 20, 2017 at 4:24 PM, Barry Smith  wrote:

>
> > On Oct 20, 2017, at 12:17 PM, Kong, Fande  wrote:
> >
> > How about to do  a global check (MPI_Allreduce)?
>
>Kill you for large number of processes.
>
> > If we do not set values at all and the matrix is already assembled, we
> just return without doing anything?
> >
> > How expensive, in the current implementation,  to call
> MatAssemblyBegin/End if there are no any stashed data? Is it so cheap that
> we can just ignore it?
>
>It requires at least one global reduction.
>

So, this global reduction is cheaper than MPI_Allreduce?  Or they are
similar.




> >
> > I am asking because we call MatAssemblyBegin/End so often in MOOSE. I
> want to make sure this is not going to bring up any performance issue.
>
>   You need to check each use in MOOSE and see WHY it is being called. If
> no reason then don't call.
>

Good point.

Fande,



>
>   Barry
>
> >
> >
> > Fande,
> >
> > On Fri, Oct 20, 2017 at 11:08 AM, Barry Smith 
> wrote:
> >
> >   One process sets a value in the matrix, the others do not. They all
> call MatAssemblyBegin(). Some processes will skip the assembly and hence
> the code will hang.
> >
> >
> > > On Oct 20, 2017, at 12:03 PM, Kong, Fande  wrote:
> > >
> > > Hi All,
> > >
> > > In Mat/VecAssemblyBegin/End, why we do not check if or not mat/vec is
> assembled. If mat/vec is already assembled, should we just return without
> doing anything?
> > >
> > > I think we have some particular reasons not to check if the matrix is
> assembled. I honestly do not know why.
> > >
> > > Thanks,
> > >
> > > Fande,
> >
> >
>
>

Re: [petsc-dev] VecView HDF5 read / write partial array

On Fri, Oct 20, 2017 at 4:32 PM, Blaise A Bourdin  wrote:

> Hi,
>
> Recent versions of exodus are implemented with netcdf4, which is based on
> hdf5.
> Instead of going through the process of rewriting exodus viewers, I would
> like to use PETSc’s HDF5 I/O operations.
> Problem is that as far as I can understand, exodus / netcdf does not use
> hdf5 time steps. Instead, it writes a single vector containing all values
> of a field at all time steps.
> (i.e. if vij denotes the i^th component of v at time step j, v is written
> as a single array v11, v21, v31, . . ., v12, v22, v32, . . ., v13, v23,  .
> . .) and so on).
> Is there a way to use VecView hdf5 to read/write only entries kn:(k+1)n-1
> in an hdf5 file?


This is exactly what HDF5 does. You write a 'hyperslab' of the array
corresponding to that time step.

   Matt


>
> Blaise
>
> --
> Department of Mathematics and Center for Computation & Technology
> Louisiana State University, Baton Rouge, LA 70803, USA
> Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~
> bourdin
>
>
>
>
>
>
>
>


-- 
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] TS Terminology


> On Oct 20, 2017, at 12:20 PM, Jed Brown  wrote:
> 
> Barry Smith  writes:
> 
>>> On Oct 20, 2017, at 11:53 AM, Jed Brown  wrote:
>>> 
>>> Barry Smith  writes:
>>> 
  The name absolutely has to be changed. But to what? And the manual page 
 is WRONG! You cannot justify that no matter how much you want to keep the 
 current confusing/inaccurate name.
>>> 
>>> Are you also adamant that SNESComputeFunction must be changed?  After
>>> all, it isn't returning the output of the function that was passed to
>>> SNESSetFunction.  If SNESComputeFunction is okay, but TSComputeIFunction
>>> is not, what is the rationale for that?
>> 
>>  In TS it is DAMN confusing. (Since you and Emil have lived with it from 
>> day one I know it is not confusing to you; but it is confusing to everyone 
>> else).
> 
> Why is it DAMN confusing in TS, but exactly the same pattern is
> "very minor" in SNES?  We need an explanation here unless we're going to
> rename all XComputeY() functions to compute exactly Y.

   Please tell us all the XComputeY() that need to be fixed so we can fix them 
all. I am ready to fix them today.

  Barry

> 
> Is it possible it has something to do with you having spent 30 years
> thinking deeply about Newton solvers and this moving the RHS vector over
> just seems like a trivial and obvious transformation?
> 
>>  In SNES it is a very minor confusion.
>> 
>>  We absolutely need to fix things that are DAMN confusing. Fixing things 
>> that are minor confusing is much less important. So it would be fine to 
>> change SNESComputeFunction() but I have no reason to be adamant about it.
>> 
>>  Barry

Re: [petsc-dev] MatAssembly and VecAssembly


> On Oct 20, 2017, at 12:17 PM, Kong, Fande  wrote:
> 
> How about to do  a global check (MPI_Allreduce)?

   Kill you for large number of processes.

> If we do not set values at all and the matrix is already assembled, we just 
> return without doing anything?  
> 
> How expensive, in the current implementation,  to call MatAssemblyBegin/End 
> if there are no any stashed data? Is it so cheap that we can just ignore it? 

   It requires at least one global reduction. 
> 
> I am asking because we call MatAssemblyBegin/End so often in MOOSE. I want to 
> make sure this is not going to bring up any performance issue. 

  You need to check each use in MOOSE and see WHY it is being called. If no 
reason then don't call. 

  Barry

> 
> 
> Fande,
> 
> On Fri, Oct 20, 2017 at 11:08 AM, Barry Smith  wrote:
> 
>   One process sets a value in the matrix, the others do not. They all call 
> MatAssemblyBegin(). Some processes will skip the assembly and hence the code 
> will hang.
> 
> 
> > On Oct 20, 2017, at 12:03 PM, Kong, Fande  wrote:
> >
> > Hi All,
> >
> > In Mat/VecAssemblyBegin/End, why we do not check if or not mat/vec is 
> > assembled. If mat/vec is already assembled, should we just return without 
> > doing anything?
> >
> > I think we have some particular reasons not to check if the matrix is 
> > assembled. I honestly do not know why.
> >
> > Thanks,
> >
> > Fande,
> 
>

Re: [petsc-dev] TS Terminology

"Zhang, Hong"  writes:

>> On Oct 20, 2017, at 3:15 PM, Jed Brown  wrote:
>> 
>> "Zhang, Hong"  writes:
>> 
>>> Another confusion which is not related to this topic is the usage of the 
>>> word "DAE".
>>> I disagree with the statement "In general, this (the general form) is a 
>>> differential algebraic equation (DAE)" on page 141 of the manual.
>>> The word "DAE" has been abused in the comments of many TS functions 
>>> (including TSComputeIFunction), where it actually should mean "DAE or ODE".
>>> PETSc uses the same interfaces for both DAE and ODE, but it is wrong to 
>>> consider ODE as a special case of DAE.
>>> They are fundamentally different from each other and should be 
>>> distinguished explicitly in the manual and the source code.
>> 
>> Is this not valid usage?
>> 
>>  "Thus, ODEs have index 0."
>
> One can think of scalar as 0d matrix. But how often do people use 0d matrix 
> for scalar?

The elegance of differential geometry would be lost if scalars were not
0-forms.

> When people say DAE, naturally it refers to DAE with index 1 or higher 
> (ODE+algebraic constraints).
> If PETSc TS developers are really DAE zealots, at least it should be 
> specified that their "DAE" includes index 0 DAE.

I'm fine with explicitly stating that ODE are index 0 DAE and having
phrases like "if your DAE is an ODE, the explicit interface can be
used".

Re: [petsc-dev] TS Terminology

2017-10-20 Thread Zhang, Hong

> On Oct 20, 2017, at 3:15 PM, Jed Brown  wrote:
> 
> "Zhang, Hong"  writes:
> 
>> Another confusion which is not related to this topic is the usage of the 
>> word "DAE".
>> I disagree with the statement "In general, this (the general form) is a 
>> differential algebraic equation (DAE)" on page 141 of the manual.
>> The word "DAE" has been abused in the comments of many TS functions 
>> (including TSComputeIFunction), where it actually should mean "DAE or ODE".
>> PETSc uses the same interfaces for both DAE and ODE, but it is wrong to 
>> consider ODE as a special case of DAE.
>> They are fundamentally different from each other and should be distinguished 
>> explicitly in the manual and the source code.
> 
> Is this not valid usage?
> 
>  "Thus, ODEs have index 0."

One can think of scalar as 0d matrix. But how often do people use 0d matrix for 
scalar?

When people say DAE, naturally it refers to DAE with index 1 or higher 
(ODE+algebraic constraints).
If PETSc TS developers are really DAE zealots, at least it should be specified 
that their "DAE" includes index 0 DAE.
Most parts of the manual distinguish well between ODE and DAE, shouldn't the 
comments in the source code be improved at least for consistency?

Hong (Mr.)

> 
>  http://www.scholarpedia.org/article/Differential-algebraic_equations

[petsc-dev] VecView HDF5 read / write partial array

2017-10-20 Thread Blaise A Bourdin

Hi,

Recent versions of exodus are implemented with netcdf4, which is based on hdf5.
Instead of going through the process of rewriting exodus viewers, I would like 
to use PETSc’s HDF5 I/O operations.
Problem is that as far as I can understand, exodus / netcdf does not use hdf5 
time steps. Instead, it writes a single vector containing all values of a field 
at all time steps.
(i.e. if vij denotes the i^th component of v at time step j, v is written as a 
single array v11, v21, v31, . . ., v12, v22, v32, . . ., v13, v23,  . . .) and 
so on).
Is there a way to use VecView hdf5 to read/write only entries kn:(k+1)n-1 in an 
hdf5 file?

Blaise

-- 
Department of Mathematics and Center for Computation & Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

Re: [petsc-dev] TS Terminology

Matthew Knepley  writes:

>> Note that TSComputeIFunction is very much like SNESComputeFunction,
>> which includes
>>
>>   if (snes->vec_rhs) {
>> ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr);
>>   }
>>
>> Why haven't you complained about that?
>>
>
> Good point. I did not notice. This came up because the initialization
> of input vectors is inconsistent between TSComputeIFunction() and
> TSComputeIFunctionLocal(). The former does not zero the output vec,
> but the later does.

The latter function doesn't exist so maybe you mean
TSComputeIFunction_DMDA with ADD_VALUES?  That's because the DMDA needs
it when using ADD_VALUES, just like SNESComputeFunction_DMDA.  When
using INSERT_VALUES, the user is responsible for setting every entry.
Is any of this different from SNES?

Re: [petsc-dev] TS Terminology

"Zhang, Hong"  writes:

> Another confusion which is not related to this topic is the usage of the word 
> "DAE".
> I disagree with the statement "In general, this (the general form) is a 
> differential algebraic equation (DAE)" on page 141 of the manual.
> The word "DAE" has been abused in the comments of many TS functions 
> (including TSComputeIFunction), where it actually should mean "DAE or ODE".
> PETSc uses the same interfaces for both DAE and ODE, but it is wrong to 
> consider ODE as a special case of DAE.
> They are fundamentally different from each other and should be distinguished 
> explicitly in the manual and the source code.

Is this not valid usage?

  "Thus, ODEs have index 0."

  http://www.scholarpedia.org/article/Differential-algebraic_equations

Re: [petsc-dev] TS Terminology

2017-10-20 Thread Zhang, Hong

On Oct 20, 2017, at 11:34 AM, Jed Brown 
> wrote:

Matthew Knepley > writes:

On Fri, Oct 20, 2017 at 11:58 AM, Emil Constantinescu 
>
wrote:

On 10/20/17 10:43 AM, Matthew Knepley wrote:

On Fri, Oct 20, 2017 at 11:22 AM, Emil Constantinescu <
emcon...@mcs.anl.gov 
> wrote:

   On 10/20/17 9:11 AM, Matthew Knepley wrote:

   On Fri, Oct 20, 2017 at 9:23 AM, Emil Constantinescu

   >>
wrote:

On 10/20/17 7:57 AM, Matthew Knepley wrote:

I am confused by some of the terminology in TS. At the
top
level, IFunction appears to mean the entire equation

F(u, u_t, x) = 0

Matt, page 141 of the manual: F(t, u, u_t) = G(t, u), and
   not zero
on the RHS side. To make the interface general we allow
   internally
for F:= F(t, u, u_t) - G(t, u) and then F=0.

   This is not "internal". Its the toplevel interface:

   https://bitbucket.org/petsc/petsc/src/63ae3ecac3af8ce782273a
76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master&
fileviewer=file-view-default#ts.c-884

   It computes F - G.

   That's what it should do in some cases. The user provides either
   ifunction or rhs funtion or both. The api to the solvers can take
   care of this stuff automatically - that's what I meant by internal.
   Different TS solvers can take different definitions of the funtions;
   e.g., imex need both, beuler can take ifuntion and/or rhs function
   but instead of writing beuler for both we choose the most general
   case (ifunction) and compose the functions accordingly. The F - G is
   transparent to the user. But somewhere the sausage needs to be made
   and I think that is the right level because that is least likely to
   change and least maintenance.

I know what it does. I looked at the code. You are missing the point here.

We cannot use the same word, IFunction, for two different things, F and
F-G. The argument that is is not user facing is complete bullshit.
The user inputs the points for ifunction, and can also call the toplevel
interface.

Matt, we do not. IFunction is F(t,u_t,u), RHS function is G(t,u). What we
solve is F=G and not F=0. Do you doubt that?

When the user specifies IFunction it is that F; when the user specifies
RHS it is that G.

Nope. We use the word IFunction, specifically in the ifunction function
pointer to mean

 F

but we use the word IFunction, specifically in TSComputeIFunction, to mean

 F - G

And, no its visible to everyone, not just "TS developers".

The "imex" flag is part of the TSComputeIFunction interface.  If that
flag is TRUE, then IFunction is exactly what the user provides.  If
FALSE, then the RHSFunction is subtracted off.  Should that be better
documented?

I think the confusion really comes from the following description in the 
comments:

TSComputeIFunction
 - Evaluates the DAE residual written in implicit form F(t,U,Udot)=0

It would be easier to understand if we say

TSComputeIFunction
 - Evaluates F or F-G, depending on the "imex" flag

since it is stated in the manual that PETSc addresses the general form 
F(t,u,udot) = G(t,u).
And normally users do not need to be aware of anything about the fully implicit 
form that some internal TS solvers work with, but developers do.

Another confusion which is not related to this topic is the usage of the word 
"DAE".
I disagree with the statement "In general, this (the general form) is a 
differential algebraic equation (DAE)" on page 141 of the manual.
The word "DAE" has been abused in the comments of many TS functions (including 
TSComputeIFunction), where it actually should mean "DAE or ODE".
PETSc uses the same interfaces for both DAE and ODE, but it is wrong to 
consider ODE as a special case of DAE.
They are fundamentally different from each other and should be distinguished 
explicitly in the manual and the source code.

Hong (Mr.)

We could make a new interface TSComputeIFunctionMaybeMinusRHSFunction()
but I don't think it's necessary.  Without any such interface, the TS
implementations would each need to reproduce a bunch of vector and
matrix wrangling.

Note that TSComputeIFunction is very much like SNESComputeFunction,
which includes

 if (snes->vec_rhs) {
   ierr =

Re: [petsc-dev] TS Terminology

On Fri, Oct 20, 2017 at 1:20 PM, Jed Brown  wrote:

> Barry Smith  writes:
>
> >> On Oct 20, 2017, at 11:53 AM, Jed Brown  wrote:
> >>
> >> Barry Smith  writes:
> >>
> >>>   The name absolutely has to be changed. But to what? And the manual
> page is WRONG! You cannot justify that no matter how much you want to keep
> the current confusing/inaccurate name.
> >>
> >> Are you also adamant that SNESComputeFunction must be changed?  After
> >> all, it isn't returning the output of the function that was passed to
> >> SNESSetFunction.  If SNESComputeFunction is okay, but TSComputeIFunction
> >> is not, what is the rationale for that?
> >
> >   In TS it is DAMN confusing. (Since you and Emil have lived with it
> from day one I know it is not confusing to you; but it is confusing to
> everyone else).
>
> Why is it DAMN confusing in TS, but exactly the same pattern is
> "very minor" in SNES?  We need an explanation here unless we're going to
> rename all XComputeY() functions to compute exactly Y.
>
> Is it possible it has something to do with you having spent 30 years
> thinking deeply about Newton solvers and this moving the RHS vector over
> just seems like a trivial and obvious transformation?


I think you are right that the SNES thing needs fixed as well. The problem
has not
come up because no one ever uses that interface. It is only there because
it makes
FAS a lot easier.

   Matt


> >   In SNES it is a very minor confusion.
> >
> >   We absolutely need to fix things that are DAMN confusing. Fixing
> things that are minor confusing is much less important. So it would be fine
> to change SNESComputeFunction() but I have no reason to be adamant about it.
> >
> >   Barry
>



-- 
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] TS Terminology

On Fri, Oct 20, 2017 at 12:34 PM, Jed Brown  wrote:

> Matthew Knepley  writes:
>
> > On Fri, Oct 20, 2017 at 11:58 AM, Emil Constantinescu <
> emcon...@mcs.anl.gov>
> > wrote:
> >
> >> On 10/20/17 10:43 AM, Matthew Knepley wrote:
> >>
> >>> On Fri, Oct 20, 2017 at 11:22 AM, Emil Constantinescu <
> >>> emcon...@mcs.anl.gov > wrote:
> >>>
> >>> On 10/20/17 9:11 AM, Matthew Knepley wrote:
> >>>
> >>> On Fri, Oct 20, 2017 at 9:23 AM, Emil Constantinescu
> >>> 
> >>> >>
> >>> wrote:
> >>>
> >>>  On 10/20/17 7:57 AM, Matthew Knepley wrote:
> >>>
> >>>  I am confused by some of the terminology in TS. At the
> >>> top
> >>>  level, IFunction appears to mean the entire equation
> >>>
> >>>  F(u, u_t, x) = 0
> >>>
> >>>
> >>>  Matt, page 141 of the manual: F(t, u, u_t) = G(t, u), and
> >>> not zero
> >>>  on the RHS side. To make the interface general we allow
> >>> internally
> >>>  for F:= F(t, u, u_t) - G(t, u) and then F=0.
> >>>
> >>>
> >>> This is not "internal". Its the toplevel interface:
> >>>
> >>> https://bitbucket.org/petsc/petsc/src/63ae3ecac3af8ce782273a
> >>> 76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master&
> >>> fileviewer=file-view-default#ts.c-884
> >>>  >>> a76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master&
> >>> fileviewer=file-view-default#ts.c-884>
> >>>
> >>> It computes F - G.
> >>>
> >>>
> >>> That's what it should do in some cases. The user provides either
> >>> ifunction or rhs funtion or both. The api to the solvers can take
> >>> care of this stuff automatically - that's what I meant by internal.
> >>> Different TS solvers can take different definitions of the
> funtions;
> >>> e.g., imex need both, beuler can take ifuntion and/or rhs function
> >>> but instead of writing beuler for both we choose the most general
> >>> case (ifunction) and compose the functions accordingly. The F - G
> is
> >>> transparent to the user. But somewhere the sausage needs to be made
> >>> and I think that is the right level because that is least likely to
> >>> change and least maintenance.
> >>>
> >>>
> >>> I know what it does. I looked at the code. You are missing the point
> here.
> >>>
> >>> We cannot use the same word, IFunction, for two different things, F and
> >>> F-G. The argument that is is not user facing is complete bullshit.
> >>> The user inputs the points for ifunction, and can also call the
> toplevel
> >>> interface.
> >>>
> >>
> >> Matt, we do not. IFunction is F(t,u_t,u), RHS function is G(t,u). What
> we
> >> solve is F=G and not F=0. Do you doubt that?
> >>
> >> When the user specifies IFunction it is that F; when the user specifies
> >> RHS it is that G.
> >>
> >
> > Nope. We use the word IFunction, specifically in the ifunction function
> > pointer to mean
> >
> >   F
> >
> > but we use the word IFunction, specifically in TSComputeIFunction, to
> mean
> >
> >   F - G
> >
> > And, no its visible to everyone, not just "TS developers".
>
> The "imex" flag is part of the TSComputeIFunction interface.  If that
> flag is TRUE, then IFunction is exactly what the user provides.  If
> FALSE, then the RHSFunction is subtracted off.  Should that be better
> documented?
>
> We could make a new interface TSComputeIFunctionMaybeMinusRHSFunction()
> but I don't think it's necessary.  Without any such interface, the TS
> implementations would each need to reproduce a bunch of vector and
> matrix wrangling.
>
> Note that TSComputeIFunction is very much like SNESComputeFunction,
> which includes
>
>   if (snes->vec_rhs) {
> ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr);
>   }
>
> Why haven't you complained about that?
>

Good point. I did not notice. This came up because the initialization of
input vectors is inconsistent
between TSComputeIFunction() and TSComputeIFunctionLocal(). The former does
not zero the
output vec, but the later does.

   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] Improving nightly builds for maint

Barry Smith  writes:

>> On Oct 18, 2017, at 3:10 PM, 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?
>
>Find memory corruption bugs that won't be found otherwise? 

So use three elements and two time steps.

>But note that the original test problems should be taking only a short 
> time so even with valgrind should not take 2 hours.

Yes.

Re: [petsc-dev] MatAssembly and VecAssembly

If you're using the same sparsity layout, you can use -vec_assembly_bts
and VecSetOption(X,VEC_SUBSET_OFF_PROC_ENTRIES) to make repeat
assemblies just send empty messages to neighbors to notify when nothing
needs to be done.  Similar with Mat SUBSET_NONZERO_PATTERN.

But MOOSE should really think through the order of operations so it only
needs one assembly to set up a system.

"Kong, Fande"  writes:

> How about to do  a global check (MPI_Allreduce)? If we do not set values at
> all and the matrix is already assembled, we just return without doing
> anything?
>
> How expensive, in the current implementation,  to call MatAssemblyBegin/End
> if there are no any stashed data? Is it so cheap that we can just ignore
> it?
>
> I am asking because we call MatAssemblyBegin/End so often in MOOSE. I want
> to make sure this is not going to bring up any performance issue.
>
>
> Fande,
>
> On Fri, Oct 20, 2017 at 11:08 AM, Barry Smith  wrote:
>
>>
>>   One process sets a value in the matrix, the others do not. They all call
>> MatAssemblyBegin(). Some processes will skip the assembly and hence the
>> code will hang.
>>
>>
>> > On Oct 20, 2017, at 12:03 PM, Kong, Fande  wrote:
>> >
>> > Hi All,
>> >
>> > In Mat/VecAssemblyBegin/End, why we do not check if or not mat/vec is
>> assembled. If mat/vec is already assembled, should we just return without
>> doing anything?
>> >
>> > I think we have some particular reasons not to check if the matrix is
>> assembled. I honestly do not know why.
>> >
>> > Thanks,
>> >
>> > Fande,
>>
>>

Re: [petsc-dev] PetscSF in Fortran

Barry Smith  writes:

>Adrian,
>
> You should not use F90Array1d *rptr as arguments in the Fortran 
> interface. You should just use a regular Fortran one dimensional array of 
> real/scalar.
> Fortran doesn't handle polymorphism in this way at all. You have to have 
> multiple f90 interfaces, one for each type and provide a C stub for real, 
> int, or whatever else you want to send.

Barry, look at the "use mpi_f08" way of calling MPI from Fortran.

>What types do you need to send? It is probably easier if we just write the 
> Fortran interfaces for you.
>
>If you can send a simple Fortran PETS  code that uses PetscSFBcastBegin() 
> and PetscSFBcastEnd() that I can use for testing then I will write them. Do a 
> bcast with int and another with real/scalar.
>
>
>Barry
>  
>
>
>> On Oct 19, 2017, at 8:59 PM, Adrian Croucher  
>> wrote:
>> 
>> hi
>> 
>> On 19/10/17 06:45, Matthew Knepley wrote:
>>> On Tue, Oct 17, 2017 at 11:35 PM, Adrian Croucher 
>>>  wrote:
>>> 
>>> 
>>> 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 think this is working now by just declaring those void* C variables as 
>> type(*) in the Fortran interface, e.g.:
>> 
>>   Interface
>>  Subroutine PetscSFBcastBegin(sf,unit,rarray,
>>  &   larray,ierr)
>>use petscisdef
>>PetscSF :: sf
>>PetscInt :: unit
>>type(*) :: rarray(:)
>>type(*) :: larray(:)
>>PetscErrorCode :: ierr
>>  End Subroutine PetscSFBcastBegin
>>   End Interface
>> 
>> The only difficulty I have left with this is in the MPI_Datatype variable. 
>> I'd forgotten that these datatypes are all different in C and Fortran as 
>> well.
>> 
>> I amended the C interface code to the following, to convert the Fortran MPI 
>> datatype (actually an integer) to a C MPI_Datatype:
>> 
>> PETSC_EXTERN void PETSC_STDCALL petscsfbcastbegin_(PetscSF *sf, MPI_Fint 
>> *unit, F90Array1d *rptr, F90Array1d *lptr , int *ierr 
>> PETSC_F90_2PTR_PROTO(rptrd) PETSC_F90_2PTR_PROTO(lptrd))
>> {
>>   PetscDataType pdtype;
>>   MPI_Datatype dtype;
>>   const void* rootdata;
>>   void* leafdata;
>> 
>>   dtype = MPI_Type_f2c(*unit);
>>   *ierr = PetscMPIDataTypeToPetscDataType(dtype, );if (*ierr) return;
>>   *ierr = F90Array1dAccess(rptr, pdtype, (void**)  
>> PETSC_F90_2PTR_PARAM(rptrd));if (*ierr) return;
>>   *ierr = F90Array1dAccess(lptr, pdtype, (void**)  
>> PETSC_F90_2PTR_PARAM(lptrd));if (*ierr) return;
>> 
>>   *ierr = PetscSFBcastBegin(*sf, dtype, rootdata, leafdata);
>> 
>> }
>> 
>> The problem is this only seems to work if I declare the datatype in the 
>> calling Fortran code to be of the appropriate C MPI datatype, e.g. MPI_INT, 
>> rather than the corresponding Fortran datatype, e.g. MPI_INTEGER (which 
>> causes PetscMPIDataTypeToPetscDataType() to fail, as something weird gets 
>> passed in for dtype).
>> 
>> I was expecting the opposite to be true. It doesn't seem right to have to 
>> use the C datatypes in Fortran code (confusing if the Fortran datatypes are 
>> used elsewhere). So I suspect I've messed something up. Anyone have any 
>> ideas?
>> 
>> - 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
>>

Re: [petsc-dev] TS Terminology

Barry Smith  writes:

>> On Oct 20, 2017, at 11:53 AM, Jed Brown  wrote:
>> 
>> Barry Smith  writes:
>> 
>>>   The name absolutely has to be changed. But to what? And the manual page 
>>> is WRONG! You cannot justify that no matter how much you want to keep the 
>>> current confusing/inaccurate name.
>> 
>> Are you also adamant that SNESComputeFunction must be changed?  After
>> all, it isn't returning the output of the function that was passed to
>> SNESSetFunction.  If SNESComputeFunction is okay, but TSComputeIFunction
>> is not, what is the rationale for that?
>
>   In TS it is DAMN confusing. (Since you and Emil have lived with it from 
> day one I know it is not confusing to you; but it is confusing to everyone 
> else).

Why is it DAMN confusing in TS, but exactly the same pattern is
"very minor" in SNES?  We need an explanation here unless we're going to
rename all XComputeY() functions to compute exactly Y.

Is it possible it has something to do with you having spent 30 years
thinking deeply about Newton solvers and this moving the RHS vector over
just seems like a trivial and obvious transformation?

>   In SNES it is a very minor confusion.
>
>   We absolutely need to fix things that are DAMN confusing. Fixing things 
> that are minor confusing is much less important. So it would be fine to 
> change SNESComputeFunction() but I have no reason to be adamant about it.
>
>   Barry

Re: [petsc-dev] MatAssembly and VecAssembly

2017-10-20 Thread Kong, Fande

How about to do  a global check (MPI_Allreduce)? If we do not set values at
all and the matrix is already assembled, we just return without doing
anything?

How expensive, in the current implementation,  to call MatAssemblyBegin/End
if there are no any stashed data? Is it so cheap that we can just ignore
it?

I am asking because we call MatAssemblyBegin/End so often in MOOSE. I want
to make sure this is not going to bring up any performance issue.

Fande,

On Fri, Oct 20, 2017 at 11:08 AM, Barry Smith  wrote:

>
>   One process sets a value in the matrix, the others do not. They all call
> MatAssemblyBegin(). Some processes will skip the assembly and hence the
> code will hang.
>
>
> > On Oct 20, 2017, at 12:03 PM, Kong, Fande  wrote:
> >
> > Hi All,
> >
> > In Mat/VecAssemblyBegin/End, why we do not check if or not mat/vec is
> assembled. If mat/vec is already assembled, should we just return without
> doing anything?
> >
> > I think we have some particular reasons not to check if the matrix is
> assembled. I honestly do not know why.
> >
> > Thanks,
> >
> > Fande,
>
>

Re: [petsc-dev] MatAssembly and VecAssembly

  One process sets a value in the matrix, the others do not. They all call 
MatAssemblyBegin(). Some processes will skip the assembly and hence the code 
will hang.

> On Oct 20, 2017, at 12:03 PM, Kong, Fande  wrote:
> 
> Hi All,
> 
> In Mat/VecAssemblyBegin/End, why we do not check if or not mat/vec is 
> assembled. If mat/vec is already assembled, should we just return without 
> doing anything?
> 
> I think we have some particular reasons not to check if the matrix is 
> assembled. I honestly do not know why. 
> 
> Thanks,
> 
> Fande,

Re: [petsc-dev] PetscSF in Fortran


   Adrian,

You should not use F90Array1d *rptr as arguments in the Fortran interface. 
You should just use a regular Fortran one dimensional array of real/scalar.
Fortran doesn't handle polymorphism in this way at all. You have to have 
multiple f90 interfaces, one for each type and provide a C stub for real, int, 
or whatever else you want to send.

   What types do you need to send? It is probably easier if we just write the 
Fortran interfaces for you.

   If you can send a simple Fortran PETS  code that uses PetscSFBcastBegin() 
and PetscSFBcastEnd() that I can use for testing then I will write them. Do a 
bcast with int and another with real/scalar.


   Barry
 


> On Oct 19, 2017, at 8:59 PM, Adrian Croucher  
> wrote:
> 
> hi
> 
> On 19/10/17 06:45, Matthew Knepley wrote:
>> On Tue, Oct 17, 2017 at 11:35 PM, Adrian Croucher 
>>  wrote:
>> 
>> 
>> 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 think this is working now by just declaring those void* C variables as 
> type(*) in the Fortran interface, e.g.:
> 
>   Interface
>  Subroutine PetscSFBcastBegin(sf,unit,rarray,
>  &   larray,ierr)
>use petscisdef
>PetscSF :: sf
>PetscInt :: unit
>type(*) :: rarray(:)
>type(*) :: larray(:)
>PetscErrorCode :: ierr
>  End Subroutine PetscSFBcastBegin
>   End Interface
> 
> The only difficulty I have left with this is in the MPI_Datatype variable. 
> I'd forgotten that these datatypes are all different in C and Fortran as well.
> 
> I amended the C interface code to the following, to convert the Fortran MPI 
> datatype (actually an integer) to a C MPI_Datatype:
> 
> PETSC_EXTERN void PETSC_STDCALL petscsfbcastbegin_(PetscSF *sf, MPI_Fint 
> *unit, F90Array1d *rptr, F90Array1d *lptr , int *ierr 
> PETSC_F90_2PTR_PROTO(rptrd) PETSC_F90_2PTR_PROTO(lptrd))
> {
>   PetscDataType pdtype;
>   MPI_Datatype dtype;
>   const void* rootdata;
>   void* leafdata;
> 
>   dtype = MPI_Type_f2c(*unit);
>   *ierr = PetscMPIDataTypeToPetscDataType(dtype, );if (*ierr) return;
>   *ierr = F90Array1dAccess(rptr, pdtype, (void**)  
> PETSC_F90_2PTR_PARAM(rptrd));if (*ierr) return;
>   *ierr = F90Array1dAccess(lptr, pdtype, (void**)  
> PETSC_F90_2PTR_PARAM(lptrd));if (*ierr) return;
> 
>   *ierr = PetscSFBcastBegin(*sf, dtype, rootdata, leafdata);
> 
> }
> 
> The problem is this only seems to work if I declare the datatype in the 
> calling Fortran code to be of the appropriate C MPI datatype, e.g. MPI_INT, 
> rather than the corresponding Fortran datatype, e.g. MPI_INTEGER (which 
> causes PetscMPIDataTypeToPetscDataType() to fail, as something weird gets 
> passed in for dtype).
> 
> I was expecting the opposite to be true. It doesn't seem right to have to use 
> the C datatypes in Fortran code (confusing if the Fortran datatypes are used 
> elsewhere). So I suspect I've messed something up. Anyone have any ideas?
> 
> - 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] MatAssembly and VecAssembly

2017-10-20 Thread Kong, Fande

Hi All,

In Mat/VecAssemblyBegin/End, why we do not check if or not mat/vec is
assembled. If mat/vec is already assembled, should we just return without
doing anything?

I think we have some particular reasons not to check if the matrix is
assembled. I honestly do not know why.

Thanks,

Fande,

Re: [petsc-dev] PetscSF in Fortran


> On Oct 19, 2017, at 10:29 PM, Jed Brown  wrote:
> 
> Adrian Croucher  writes:
> 
>> hi
>> 
>> 
>> The problem is this only seems to work if I declare the datatype in the 
>> calling Fortran code to be of the appropriate C MPI datatype, e.g. 
>> MPI_INT, rather than the corresponding Fortran datatype, e.g. 
>> MPI_INTEGER (which causes PetscMPIDataTypeToPetscDataType() to fail, as 
>> something weird gets passed in for dtype).
>> 
>> I was expecting the opposite to be true. It doesn't seem right to have 
>> to use the C datatypes in Fortran code (confusing if the Fortran 
>> datatypes are used elsewhere). So I suspect I've messed something up. 
>> Anyone have any ideas?
> 
> I think PetscMPIDataTypeToPetscDataType should be extended to handle the
> Fortran types.
> 
> BTW, this code looks very wrong to me:
> 
>  if (mtype == MPIU_INT) *ptype = PETSC_INT;
>  else if (mtype == MPI_INT) *ptype = PETSC_INT;

   Yes definitely wrong. 

> 
> [...]
> 
>  else if (mtype == MPI_LONG)*ptype = PETSC_LONG;
> 
> 
> If PetscInt is 64-bit then MPI_INT (which corresponds to int) would be
> mapped to a 64-bit int.  Those PETSc datatypes either need to be
> abandoned (I think MPI types would do) or we need a terminology to
> distinguish PetscInt from int.

  From a manual page

-
/*E
PetscDataType - Used for handling different basic data types.

   Level: beginner

   Developer comment: It would be nice if we could always just use MPI 
Datatypes, why can we not?
--

   When MPI first came out I struggled to use just MPI datatypes and not 
introduce a parallel set of PETSc ones; I was unsuccessful but that may have 
just have come from my ignorance and not because it was not possible.

   Regardless of whether we can totally eliminate the use of PETSc datatypes we 
should try to limit them and NOT use them when not needed.


I have made two pull requests 

1) fix the bug 
https://bitbucket.org/petsc/petsc/pull-requests/782/mpi_int-mpi-datatype-with-64-bit-int/diff

2) removed unneeded use of PetscDataType 
https://bitbucket.org/petsc/petsc/pull-requests/783/remove-use-of-petscdatatype-from-vtk-code/diff

Note this does not resolve Adrian's problem. I will respond to Adrian in 
another email.

  Barry

Re: [petsc-dev] TS Terminology

> On Oct 20, 2017, at 11:56 AM, Emil Constantinescu  
> wrote:
> 
> 
> 
> On 10/20/17 11:39 AM, Barry Smith wrote:
>>You guys are arguing about irrelevancies. Focus on the NAME, that is what 
>> Matt is complaining about.
> 
> How about we change the manual? Isn't that easier?

   The manual of course needs to be fixed, but that is not enough. The function 
name is fundamentally misleading and thus needs to be fixed.

  Barry

> 
> Emil

Re: [petsc-dev] TS Terminology

> On Oct 20, 2017, at 11:53 AM, Jed Brown  wrote:
> 
> Barry Smith  writes:
> 
>>   The name absolutely has to be changed. But to what? And the manual page is 
>> WRONG! You cannot justify that no matter how much you want to keep the 
>> current confusing/inaccurate name.
> 
> Are you also adamant that SNESComputeFunction must be changed?  After
> all, it isn't returning the output of the function that was passed to
> SNESSetFunction.  If SNESComputeFunction is okay, but TSComputeIFunction
> is not, what is the rationale for that?

  In TS it is DAMN confusing. (Since you and Emil have lived with it from 
day one I know it is not confusing to you; but it is confusing to everyone 
else).

  In SNES it is a very minor confusion.

  We absolutely need to fix things that are DAMN confusing. Fixing things 
that are minor confusing is much less important. So it would be fine to change 
SNESComputeFunction() but I have no reason to be adamant about it.

  Barry

Re: [petsc-dev] TS Terminology




On 10/20/17 11:39 AM, Barry Smith wrote:

You guys are arguing about irrelevancies. Focus on the NAME, that is what 
Matt is complaining about.


How about we change the manual? Isn't that easier?

Emil

Re: [petsc-dev] TS Terminology

Barry Smith  writes:

>The name absolutely has to be changed. But to what? And the manual page is 
> WRONG! You cannot justify that no matter how much you want to keep the 
> current confusing/inaccurate name.

Are you also adamant that SNESComputeFunction must be changed?  After
all, it isn't returning the output of the function that was passed to
SNESSetFunction.  If SNESComputeFunction is okay, but TSComputeIFunction
is not, what is the rationale for that?

Re: [petsc-dev] TS Terminology


   You guys are arguing about irrelevancies. Focus on the NAME, that is what 
Matt is complaining about.


> On Oct 20, 2017, at 11:37 AM, Jed Brown  wrote:
> 
> Emil Constantinescu  writes:
> 
>> Matt, that depends, if TS method is imex, 
> 
> Nothing to do with "TS method".  The "imex" parameter is an argument to
> the function TSComputeIFunction.  SNESComputeFunction's behavior depends
> on whether snes->vec_rhs exists, but TSComputeIFunction's behavior
> depends only on the arguments that are passed.
> 
>> then it computes just F, not F-G so your argument is not correct. If
>> the method can do only implicit it computes F and subtracts G *if
>> defined*. If the TS method can only do explicit and you define F then
>> it fails.
>> 
>> Again, this has to do with the TS methods and PETSc doing the work for 
>> you of packing the functions in different ways.
>

Re: [petsc-dev] TS Terminology


> On Oct 20, 2017, at 11:35 AM, Emil Constantinescu  
> wrote:
> 
> To me IFunction really means IFunction obtained from implicit part  (- 
> potentially explicit part depending on the solver). It is not just the 
> implicit part, it is the implicit function (not part). The default is just 
> the part thing. To me this makes it simple, less interfaces and it is 
> transparent to the user. This will change the interfaces in all solvers, will 
> fix little and make up really long names along the way. I mean this may 
> confuse users by implying that there are multiple parts whereas most users 
> have just one.
> 
> My point is the following: how does this affect how Matt implements things? 
> If you two find this is really worth it, I have nothing against it. To me 
> it's just a name.

   The name absolutely has to be changed. But to what? And the manual page is 
WRONG! You cannot justify that no matter how much you want to keep the current 
confusing/inaccurate name.

   Barry


> 
> Emil
> 
> 
> On 10/20/17 11:22 AM, Barry Smith wrote:
>>   Emil, see my email just sent. We need to rename this function (and its 
>> Jacobian partner).
>>> On Oct 20, 2017, at 11:20 AM, Emil Constantinescu  
>>> wrote:
>>> 
>>> 
>>> 
>>> Matt, that depends, if TS method is imex, then it computes just F, not F-G 
>>> so your argument is not correct. If the method can do only implicit it 
>>> computes F and subtracts G *if defined*. If the TS method can only do 
>>> explicit and you define F then it fails.
>>> 
>>> Again, this has to do with the TS methods and PETSc doing the work for you 
>>> of packing the functions in different ways.
>>> 
>>> Emil
>>> 
   Matt
Now internally, because different solvers have different needs the
IFunction ... presented to the TS solver may look differently. This
is a design choice - if you are not a TS developer it should not
affect you.
This is a design decision: if implemented at this level, we avoid
having every TS method be aware of the upper level functions.
Emil
 -- 
 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] TS Terminology

Emil Constantinescu  writes:

> Matt, that depends, if TS method is imex, 

Nothing to do with "TS method".  The "imex" parameter is an argument to
the function TSComputeIFunction.  SNESComputeFunction's behavior depends
on whether snes->vec_rhs exists, but TSComputeIFunction's behavior
depends only on the arguments that are passed.

> then it computes just F, not F-G so your argument is not correct. If
> the method can do only implicit it computes F and subtracts G *if
> defined*. If the TS method can only do explicit and you define F then
> it fails.
>
> Again, this has to do with the TS methods and PETSc doing the work for 
> you of packing the functions in different ways.

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


> On Oct 18, 2017, at 3:10 PM, 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?

   Find memory corruption bugs that won't be found otherwise? 

   But note that the original test problems should be taking only a short time 
so even with valgrind should not take 2 hours.

   Barry

> 
> 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] TS Terminology

To me IFunction really means IFunction obtained from implicit part  (- 
potentially explicit part depending on the solver). It is not just the 
implicit part, it is the implicit function (not part). The default is 
just the part thing. To me this makes it simple, less interfaces and it 
is transparent to the user. This will change the interfaces in all 
solvers, will fix little and make up really long names along the way. I 
mean this may confuse users by implying that there are multiple parts 
whereas most users have just one.


My point is the following: how does this affect how Matt implements 
things? If you two find this is really worth it, I have nothing against 
it. To me it's just a name.


Emil


On 10/20/17 11:22 AM, Barry Smith wrote:


   Emil, see my email just sent. We need to rename this function (and its 
Jacobian partner).


On Oct 20, 2017, at 11:20 AM, Emil Constantinescu  wrote:



Matt, that depends, if TS method is imex, then it computes just F, not F-G so 
your argument is not correct. If the method can do only implicit it computes F 
and subtracts G *if defined*. If the TS method can only do explicit and you 
define F then it fails.

Again, this has to do with the TS methods and PETSc doing the work for you of 
packing the functions in different ways.

Emil


   Matt
Now internally, because different solvers have different needs the
IFunction ... presented to the TS solver may look differently. This
is a design choice - if you are not a TS developer it should not
affect you.
This is a design decision: if implemented at this level, we avoid
having every TS method be aware of the upper level functions.
Emil
--
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] TS Terminology

Matthew Knepley  writes:

> On Fri, Oct 20, 2017 at 11:58 AM, Emil Constantinescu 
> wrote:
>
>> On 10/20/17 10:43 AM, Matthew Knepley wrote:
>>
>>> On Fri, Oct 20, 2017 at 11:22 AM, Emil Constantinescu <
>>> emcon...@mcs.anl.gov > wrote:
>>>
>>> On 10/20/17 9:11 AM, Matthew Knepley wrote:
>>>
>>> On Fri, Oct 20, 2017 at 9:23 AM, Emil Constantinescu
>>> 
>>> >>
>>> wrote:
>>>
>>>  On 10/20/17 7:57 AM, Matthew Knepley wrote:
>>>
>>>  I am confused by some of the terminology in TS. At the
>>> top
>>>  level, IFunction appears to mean the entire equation
>>>
>>>  F(u, u_t, x) = 0
>>>
>>>
>>>  Matt, page 141 of the manual: F(t, u, u_t) = G(t, u), and
>>> not zero
>>>  on the RHS side. To make the interface general we allow
>>> internally
>>>  for F:= F(t, u, u_t) - G(t, u) and then F=0.
>>>
>>>
>>> This is not "internal". Its the toplevel interface:
>>>
>>> https://bitbucket.org/petsc/petsc/src/63ae3ecac3af8ce782273a
>>> 76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master&
>>> fileviewer=file-view-default#ts.c-884
>>> >> a76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master&
>>> fileviewer=file-view-default#ts.c-884>
>>>
>>> It computes F - G.
>>>
>>>
>>> That's what it should do in some cases. The user provides either
>>> ifunction or rhs funtion or both. The api to the solvers can take
>>> care of this stuff automatically - that's what I meant by internal.
>>> Different TS solvers can take different definitions of the funtions;
>>> e.g., imex need both, beuler can take ifuntion and/or rhs function
>>> but instead of writing beuler for both we choose the most general
>>> case (ifunction) and compose the functions accordingly. The F - G is
>>> transparent to the user. But somewhere the sausage needs to be made
>>> and I think that is the right level because that is least likely to
>>> change and least maintenance.
>>>
>>>
>>> I know what it does. I looked at the code. You are missing the point here.
>>>
>>> We cannot use the same word, IFunction, for two different things, F and
>>> F-G. The argument that is is not user facing is complete bullshit.
>>> The user inputs the points for ifunction, and can also call the toplevel
>>> interface.
>>>
>>
>> Matt, we do not. IFunction is F(t,u_t,u), RHS function is G(t,u). What we
>> solve is F=G and not F=0. Do you doubt that?
>>
>> When the user specifies IFunction it is that F; when the user specifies
>> RHS it is that G.
>>
>
> Nope. We use the word IFunction, specifically in the ifunction function
> pointer to mean
>
>   F
>
> but we use the word IFunction, specifically in TSComputeIFunction, to mean
>
>   F - G
>
> And, no its visible to everyone, not just "TS developers".

The "imex" flag is part of the TSComputeIFunction interface.  If that
flag is TRUE, then IFunction is exactly what the user provides.  If
FALSE, then the RHSFunction is subtracted off.  Should that be better
documented?

We could make a new interface TSComputeIFunctionMaybeMinusRHSFunction()
but I don't think it's necessary.  Without any such interface, the TS
implementations would each need to reproduce a bunch of vector and
matrix wrangling.

Note that TSComputeIFunction is very much like SNESComputeFunction,
which includes

  if (snes->vec_rhs) {
ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr);
  }

Why haven't you complained about that?

Re: [petsc-dev] TS Terminology


  Emil, see my email just sent. We need to rename this function (and its 
Jacobian partner).

> On Oct 20, 2017, at 11:20 AM, Emil Constantinescu  
> wrote:
> 
> 
> 
> Matt, that depends, if TS method is imex, then it computes just F, not F-G so 
> your argument is not correct. If the method can do only implicit it computes 
> F and subtracts G *if defined*. If the TS method can only do explicit and you 
> define F then it fails.
> 
> Again, this has to do with the TS methods and PETSc doing the work for you of 
> packing the functions in different ways.
> 
> Emil
> 
>>   Matt
>>Now internally, because different solvers have different needs the
>>IFunction ... presented to the TS solver may look differently. This
>>is a design choice - if you are not a TS developer it should not
>>affect you.
>>This is a design decision: if implemented at this level, we avoid
>>having every TS method be aware of the upper level functions.
>>Emil
>> -- 
>> 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] TS Terminology


  Emil,

Matt is correct. He is complaining about the NAME of the function 
TSComputeIFunction; and he should also complain about its manual page that 
states

TSComputeIFunction - Evaluates the DAE residual written in implicit form 
F(t,U,Udot)=0

When imex flag is on the function computes F(t,U,Udot) but when imex flag is 
off it computes F(t,U,Udot) - G(t,U)

I think he is requesting this function name be changed (and the docs be fixed). 
Presumably similarly for the Jacobian also.

A better name would be TSComputeImplicitPartofFunction or even better 
TSComputeImplicitlyHandledPartofFunction.  The name TSComputeIFunction is 
fundamentally confusing (it is not to you because in your head you are 
translating it mentally to TSComputeImplicitlyHandledPartofFunction()

  Now does anyone have better names for this and are there other similar name 
problems that need to be fixed? Perhaps a bunch of other names also could 
benefit from renaming?

  Barry





> On Oct 20, 2017, at 11:06 AM, Matthew Knepley  wrote:
> 
> On Fri, Oct 20, 2017 at 11:58 AM, Emil Constantinescu  
> wrote:
> On 10/20/17 10:43 AM, Matthew Knepley wrote:
> On Fri, Oct 20, 2017 at 11:22 AM, Emil Constantinescu  > wrote:
> 
> On 10/20/17 9:11 AM, Matthew Knepley wrote:
> 
> On Fri, Oct 20, 2017 at 9:23 AM, Emil Constantinescu
> 
> >> wrote:
> 
>  On 10/20/17 7:57 AM, Matthew Knepley wrote:
> 
>  I am confused by some of the terminology in TS. At the top
>  level, IFunction appears to mean the entire equation
> 
>  F(u, u_t, x) = 0
> 
> 
>  Matt, page 141 of the manual: F(t, u, u_t) = G(t, u), and
> not zero
>  on the RHS side. To make the interface general we allow
> internally
>  for F:= F(t, u, u_t) - G(t, u) and then F=0.
> 
> 
> This is not "internal". Its the toplevel interface:
> 
> 
> https://bitbucket.org/petsc/petsc/src/63ae3ecac3af8ce782273a76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master=file-view-default#ts.c-884
> 
> 
> 
> It computes F - G.
> 
> 
> That's what it should do in some cases. The user provides either
> ifunction or rhs funtion or both. The api to the solvers can take
> care of this stuff automatically - that's what I meant by internal.
> Different TS solvers can take different definitions of the funtions;
> e.g., imex need both, beuler can take ifuntion and/or rhs function
> but instead of writing beuler for both we choose the most general
> case (ifunction) and compose the functions accordingly. The F - G is
> transparent to the user. But somewhere the sausage needs to be made
> and I think that is the right level because that is least likely to
> change and least maintenance.
> 
> 
> I know what it does. I looked at the code. You are missing the point here.
> 
> We cannot use the same word, IFunction, for two different things, F and F-G. 
> The argument that is is not user facing is complete bullshit.
> The user inputs the points for ifunction, and can also call the toplevel 
> interface.
> 
> Matt, we do not. IFunction is F(t,u_t,u), RHS function is G(t,u). What we 
> solve is F=G and not F=0. Do you doubt that?
> 
> When the user specifies IFunction it is that F; when the user specifies RHS 
> it is that G.
> 
> Nope. We use the word IFunction, specifically in the ifunction function 
> pointer to mean
> 
>   F
> 
> but we use the word IFunction, specifically in TSComputeIFunction, to mean
> 
>   F - G
> 
> And, no its visible to everyone, not just "TS developers".
> 
>   Matt
>  
> Now internally, because different solvers have different needs the IFunction 
> ... presented to the TS solver may look differently. This is a design choice 
> - if you are not a TS developer it should not affect you.
> 
> This is a design decision: if implemented at this level, we avoid having 
> every TS method be aware of the upper level functions.
> 
> Emil
> 
> 
> 
> 
> 
> 
> -- 
> 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] TS Terminology

On 10/20/17 11:06 AM, Matthew Knepley wrote:
On Fri, Oct 20, 2017 at 11:58 AM, Emil Constantinescu 
> wrote:

On 10/20/17 10:43 AM, Matthew Knepley wrote:

On Fri, Oct 20, 2017 at 11:22 AM, Emil Constantinescu

>> wrote:

     On 10/20/17 9:11 AM, Matthew Knepley wrote:

         On Fri, Oct 20, 2017 at 9:23 AM, Emil Constantinescu

>

>

         It computes F - G.

     That's what it should do in some cases. The user provides
either
     ifunction or rhs funtion or both. The api to the solvers
can take
     care of this stuff automatically - that's what I meant by
internal.
     Different TS solvers can take different definitions of the
funtions;
     e.g., imex need both, beuler can take ifuntion and/or rhs
function
     but instead of writing beuler for both we choose the most
general
     case (ifunction) and compose the functions accordingly. The
F - G is
     transparent to the user. But somewhere the sausage needs to
be made
     and I think that is the right level because that is least
likely to
     change and least maintenance.

I know what it does. I looked at the code. You are missing the
point here.

We cannot use the same word, IFunction, for two different
things, F and F-G. The argument that is is not user facing is
complete bullshit.
The user inputs the points for ifunction, and can also call the
toplevel interface.

Matt, we do not. IFunction is F(t,u_t,u), RHS function is G(t,u).
What we solve is F=G and not F=0. Do you doubt that?

When the user specifies IFunction it is that F; when the user
specifies RHS it is that G.

Nope. We use the word IFunction, specifically in the ifunction function 
pointer to mean

   F

but we use the word IFunction, specifically in TSComputeIFunction, to mean

   F - G

And, no its visible to everyone, not just "TS developers".

Matt, that depends, if TS method is imex, then it computes just F, not 
F-G so your argument is not correct. If the method can do only implicit 
it computes F and subtracts G *if defined*. If the TS method can only do 
explicit and you define F then it fails.

Again, this has to do with the TS methods and PETSc doing the work for 
you of packing the functions in different ways.

Emil

   Matt

Now internally, because different solvers have different needs the
IFunction ... presented to the TS solver may look differently. This
is a design choice - if you are not a TS developer it should not
affect you.

This is a design decision: if implemented at this level, we avoid
having every TS method be aware of the upper level functions.

Emil

--
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] TS Terminology

On Fri, Oct 20, 2017 at 11:58 AM, Emil Constantinescu 
wrote:

> On 10/20/17 10:43 AM, Matthew Knepley wrote:
>
>> On Fri, Oct 20, 2017 at 11:22 AM, Emil Constantinescu <
>> emcon...@mcs.anl.gov > wrote:
>>
>> On 10/20/17 9:11 AM, Matthew Knepley wrote:
>>
>> On Fri, Oct 20, 2017 at 9:23 AM, Emil Constantinescu
>> 
>> >>
>> wrote:
>>
>>  On 10/20/17 7:57 AM, Matthew Knepley wrote:
>>
>>  I am confused by some of the terminology in TS. At the
>> top
>>  level, IFunction appears to mean the entire equation
>>
>>  F(u, u_t, x) = 0
>>
>>
>>  Matt, page 141 of the manual: F(t, u, u_t) = G(t, u), and
>> not zero
>>  on the RHS side. To make the interface general we allow
>> internally
>>  for F:= F(t, u, u_t) - G(t, u) and then F=0.
>>
>>
>> This is not "internal". Its the toplevel interface:
>>
>> https://bitbucket.org/petsc/petsc/src/63ae3ecac3af8ce782273a
>> 76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master&
>> fileviewer=file-view-default#ts.c-884
>> > a76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master&
>> fileviewer=file-view-default#ts.c-884>
>>
>> It computes F - G.
>>
>>
>> That's what it should do in some cases. The user provides either
>> ifunction or rhs funtion or both. The api to the solvers can take
>> care of this stuff automatically - that's what I meant by internal.
>> Different TS solvers can take different definitions of the funtions;
>> e.g., imex need both, beuler can take ifuntion and/or rhs function
>> but instead of writing beuler for both we choose the most general
>> case (ifunction) and compose the functions accordingly. The F - G is
>> transparent to the user. But somewhere the sausage needs to be made
>> and I think that is the right level because that is least likely to
>> change and least maintenance.
>>
>>
>> I know what it does. I looked at the code. You are missing the point here.
>>
>> We cannot use the same word, IFunction, for two different things, F and
>> F-G. The argument that is is not user facing is complete bullshit.
>> The user inputs the points for ifunction, and can also call the toplevel
>> interface.
>>
>
> Matt, we do not. IFunction is F(t,u_t,u), RHS function is G(t,u). What we
> solve is F=G and not F=0. Do you doubt that?
>
> When the user specifies IFunction it is that F; when the user specifies
> RHS it is that G.
>

Nope. We use the word IFunction, specifically in the ifunction function
pointer to mean

  F

but we use the word IFunction, specifically in TSComputeIFunction, to mean

  F - G

And, no its visible to everyone, not just "TS developers".

  Matt


> Now internally, because different solvers have different needs the
> IFunction ... presented to the TS solver may look differently. This is a
> design choice - if you are not a TS developer it should not affect you.
>
> This is a design decision: if implemented at this level, we avoid having
> every TS method be aware of the upper level functions.
>
> Emil
>
>
>
>


-- 
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] TS Terminology

On 10/20/17 10:43 AM, Matthew Knepley wrote:
On Fri, Oct 20, 2017 at 11:22 AM, Emil Constantinescu
> wrote:

On 10/20/17 9:11 AM, Matthew Knepley wrote:

On Fri, Oct 20, 2017 at 9:23 AM, Emil Constantinescu

>> wrote:

On 10/20/17 7:57 AM, Matthew Knepley wrote:

I am confused by some of the terminology in TS. At the top
level, IFunction appears to mean the entire equation

F(u, u_t, x) = 0

Matt, page 141 of the manual: F(t, u, u_t) = G(t, u), and
not zero
on the RHS side. To make the interface general we allow
internally
for F:= F(t, u, u_t) - G(t, u) and then F=0.

This is not "internal". Its the toplevel interface:

https://bitbucket.org/petsc/petsc/src/63ae3ecac3af8ce782273a76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master=file-view-default#ts.c-884

It computes F - G.

That's what it should do in some cases. The user provides either
ifunction or rhs funtion or both. The api to the solvers can take
care of this stuff automatically - that's what I meant by internal.
Different TS solvers can take different definitions of the funtions;
e.g., imex need both, beuler can take ifuntion and/or rhs function
but instead of writing beuler for both we choose the most general
case (ifunction) and compose the functions accordingly. The F - G is
transparent to the user. But somewhere the sausage needs to be made
and I think that is the right level because that is least likely to
change and least maintenance.

I know what it does. I looked at the code. You are missing the point here.

We cannot use the same word, IFunction, for two different things, F and
F-G. The argument that is is not user facing is complete bullshit.
The user inputs the points for ifunction, and can also call the toplevel
interface.

Matt, we do not. IFunction is F(t,u_t,u), RHS function is G(t,u). What
we solve is F=G and not F=0. Do you doubt that?

When the user specifies IFunction it is that F; when the user specifies
RHS it is that G.

Now internally, because different solvers have different needs the
IFunction ... presented to the TS solver may look differently. This is a
design choice - if you are not a TS developer it should not affect you.

This is a design decision: if implemented at this level, we avoid having
every TS method be aware of the upper level functions.

Emil

Re: [petsc-dev] TS Terminology

On Fri, Oct 20, 2017 at 11:22 AM, Emil Constantinescu 
wrote:

> On 10/20/17 9:11 AM, Matthew Knepley wrote:
>
>> On Fri, Oct 20, 2017 at 9:23 AM, Emil Constantinescu <
>> emcon...@mcs.anl.gov > wrote:
>>
>> On 10/20/17 7:57 AM, Matthew Knepley wrote:
>>
>> I am confused by some of the terminology in TS. At the top
>> level, IFunction appears to mean the entire equation
>>
>> F(u, u_t, x) = 0
>>
>>
>> Matt, page 141 of the manual: F(t, u, u_t) = G(t, u), and not zero
>> on the RHS side. To make the interface general we allow internally
>> for F:= F(t, u, u_t) - G(t, u) and then F=0.
>>
>>
>> This is not "internal". Its the toplevel interface:
>>
>> https://bitbucket.org/petsc/petsc/src/63ae3ecac3af8ce782273a
>> 76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master&
>> fileviewer=file-view-default#ts.c-884
>>
>> It computes F - G.
>>
>
> That's what it should do in some cases. The user provides either ifunction
> or rhs funtion or both. The api to the solvers can take care of this stuff
> automatically - that's what I meant by internal. Different TS solvers can
> take different definitions of the funtions; e.g., imex need both, beuler
> can take ifuntion and/or rhs function but instead of writing beuler for
> both we choose the most general case (ifunction) and compose the functions
> accordingly. The F - G is transparent to the user. But somewhere the
> sausage needs to be made and I think that is the right level because that
> is least likely to change and least maintenance.
>

I know what it does. I looked at the code. You are missing the point here.

We cannot use the same word, IFunction, for two different things, F and
F-G. The argument that is is not user facing is complete bullshit.
The user inputs the points for ifunction, and can also call the toplevel
interface.

   Matt


> Emil
>
> Matt
>>
>> Emil
>>
>> PS: Manual: PETSc 3.8 September 26, 2017
>>
>> However, this appears to mean something different than ifunction
>> at the function pointer level. Inside TSCompiteIFunction(), it
>> uses both ifunction and rhsfunction. This makes it hard to
>> understand how composition works. TSComputeRHS() is called
>> inside TSComputeIFunction(), so if we want to reuse vectors it
>> should not initialize the vector, but it seems like
>> TSComputeIFunction() should initialize the vector since
>> TSComputeIFunctionLocal() does.
>>
>> What guarantees about initialization should we have?
>>
>>  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/
>> 
>> > >>
>>
>>
>>
>>
>> --
>> 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/ 
>>
>


-- 
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] TS Terminology

On 10/20/17 9:11 AM, Matthew Knepley wrote:
On Fri, Oct 20, 2017 at 9:23 AM, Emil Constantinescu
> wrote:

On 10/20/17 7:57 AM, Matthew Knepley wrote:

I am confused by some of the terminology in TS. At the top
level, IFunction appears to mean the entire equation

F(u, u_t, x) = 0

Matt, page 141 of the manual: F(t, u, u_t) = G(t, u), and not zero
on the RHS side. To make the interface general we allow internally
for F:= F(t, u, u_t) - G(t, u) and then F=0.

This is not "internal". Its the toplevel interface:

https://bitbucket.org/petsc/petsc/src/63ae3ecac3af8ce782273a76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master=file-view-default#ts.c-884

It computes F - G.

That's what it should do in some cases. The user provides either
ifunction or rhs funtion or both. The api to the solvers can take care
of this stuff automatically - that's what I meant by internal. Different
TS solvers can take different definitions of the funtions; e.g., imex
need both, beuler can take ifuntion and/or rhs function but instead of
writing beuler for both we choose the most general case (ifunction) and
compose the functions accordingly. The F - G is transparent to the user.
But somewhere the sausage needs to be made and I think that is the right
level because that is least likely to change and least maintenance.

Emil

Matt

Emil

PS: Manual: PETSc 3.8 September 26, 2017

However, this appears to mean something different than ifunction
at the function pointer level. Inside TSCompiteIFunction(), it
uses both ifunction and rhsfunction. This makes it hard to
understand how composition works. TSComputeRHS() is called
inside TSComputeIFunction(), so if we want to reuse vectors it
should not initialize the vector, but it seems like
TSComputeIFunction() should initialize the vector since
TSComputeIFunctionLocal() does.

What guarantees about initialization should we have?

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/

--
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] TS Terminology

On Fri, Oct 20, 2017 at 9:23 AM, Emil Constantinescu 
wrote:

> On 10/20/17 7:57 AM, Matthew Knepley wrote:
>
>> I am confused by some of the terminology in TS. At the top level,
>> IFunction appears to mean the entire equation
>>
>>F(u, u_t, x) = 0
>>
>
> Matt, page 141 of the manual: F(t, u, u_t) = G(t, u), and not zero on the
> RHS side. To make the interface general we allow internally for F:= F(t, u,
> u_t) - G(t, u) and then F=0.
>

This is not "internal". Its the toplevel interface:


https://bitbucket.org/petsc/petsc/src/63ae3ecac3af8ce782273a76ad4152cddc2fd80a/src/ts/interface/ts.c?at=master=file-view-default#ts.c-884

It computes F - G.

   Matt


> Emil
>
> PS: Manual: PETSc 3.8 September 26, 2017
>
> However, this appears to mean something different than ifunction at the
>> function pointer level. Inside TSCompiteIFunction(), it uses both ifunction
>> and rhsfunction. This makes it hard to understand how composition works.
>> TSComputeRHS() is called inside TSComputeIFunction(), so if we want to
>> reuse vectors it should not initialize the vector, but it seems like
>> TSComputeIFunction() should initialize the vector since
>> TSComputeIFunctionLocal() does.
>>
>> What guarantees about initialization should we have?
>>
>> 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/ 
>>
>


-- 
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] TS Terminology


On 10/20/17 7:57 AM, Matthew Knepley wrote:
I am confused by some of the terminology in TS. At the top level, 
IFunction appears to mean the entire equation


   F(u, u_t, x) = 0


Matt, page 141 of the manual: F(t, u, u_t) = G(t, u), and not zero on 
the RHS side. To make the interface general we allow internally for F:= 
F(t, u, u_t) - G(t, u) and then F=0.


Emil

PS: Manual: PETSc 3.8 September 26, 2017

However, this appears to mean something different than ifunction at the 
function pointer level. Inside TSCompiteIFunction(), it uses both 
ifunction and rhsfunction. This makes it hard to understand how 
composition works. TSComputeRHS() is called inside TSComputeIFunction(), 
so if we want to reuse vectors it should not initialize the vector, but 
it seems like TSComputeIFunction() should initialize the vector since 
TSComputeIFunctionLocal() does.


What guarantees about initialization should we have?

    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/

[petsc-dev] TS Terminology