[petsc-users] TS tutorial ex11 in Fortran
Jed Brown
jed at jedbrown.org
Mon Dec 28 08:30:23 CST 2020
Thibault Bridel-Bertomeu <thibault.bridelbertomeu at gmail.com> writes:
> Good morning everyone,
>
> Thank you Barry for the answer, it works now !
>
> I am facing (yet) another situation: the TSAdaptRegister function.
> In the MR on gitlab, Jed mentioned that sometimes, when function pointers
> are not stored in PETSc objects, one can use stack memory to pass that
> pointer from fortran to C.
The issue with stack memory is that when it returns, that memory is invalid. You can't use it in this instance.
I think you're going to have problems implementing a TSAdaptCreate_XYZ in Fortran (because the body of that function will need to access private struct members; see below).
I would implement what you need in C and you can call out to Fortran if you want from inside TSAdaptChoose_YourMethod().
PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt)
{
TSAdapt_DSP *dsp;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscNewLog(adapt,&dsp);CHKERRQ(ierr);
adapt->reject_safety = 1.0; /* unused */
adapt->data = (void*)dsp;
adapt->ops->choose = TSAdaptChoose_DSP;
adapt->ops->setfromoptions = TSAdaptSetFromOptions_DSP;
adapt->ops->destroy = TSAdaptDestroy_DSP;
adapt->ops->view = TSAdaptView_DSP;
ierr = PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetFilter_C",TSAdaptDSPSetFilter_DSP);CHKERRQ(ierr);
ierr = PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetPID_C",TSAdaptDSPSetPID_DSP);CHKERRQ(ierr);
ierr = TSAdaptDSPSetFilter_DSP(adapt,"PI42");CHKERRQ(ierr);
ierr = TSAdaptRestart_DSP(adapt);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
> Can anyone develop that idea ? Because for TSAdaptRegister, i guess the
> wrapper would start like :
>
> PETSC_EXTERN void tsadaptregister_(char *sname,
> void (*func)(TSAdapt*,PetscErrorCode*),
> PetscErrorCode *ierr,
> PETSC_FORTRAN_CHARLEN_T snamelen)
>
> but then the C TSAdaptRegister function takes a PetscErrorCode
> (*func)(TSAdapt) function pointer as argument ... I cannot use any
> FORTRAN_CALLBACK here since I do not have any object to hook it to, and I
> could not find a similar situation among the pre-existing wrappers. Does
> anyone have an idea on how to proceed ?
>
> Thanks !!
>
> Thibault
>
> Le mar. 22 déc. 2020 à 21:20, Barry Smith <bsmith at petsc.dev> a écrit :
>
>>
>> PetscObjectUseFortranCallback((PetscDS)ctx,
>>
>>
>> *ierr = PetscObjectSetFortranCallback((PetscObject)*prob
>>
>>
>> It looks like the problem is that these user provided functions do not
>> take a PetscDS directly as an argument so the Fortran callback information
>> cannot be obtained from them.
>>
>> The manual page for PetscDSAddBoundary() says
>>
>> - ctx - An optional user context for bcFunc
>>
>> but then when it lists the calling sequence for bcFunc it does not list
>> the ctx as an argument, so either the manual page or code is wrong.
>>
>> It looks like you make the ctx be the PetscDS prob argument when you
>> call PetscDSAddBoundary
>>
>> In principle this sounds like it might work. I think you need to track
>> through the debugger to see if the ctx passed to ourbocofunc() is
>> actually the PetscDS prob variable and if not why it is not.
>>
>> Barry
>>
>>
>> On Dec 22, 2020, at 5:49 AM, Thibault Bridel-Bertomeu <
>> thibault.bridelbertomeu at gmail.com> wrote:
>>
>> Dear all,
>>
>> I have hit two snags while implementing the missing wrappers necessary to
>> transcribe ex11 to Fortran.
>>
>> First is about the PetscDSAddBoundary wrapper, that I have done so :
>>
>> static PetscErrorCode ourbocofunc(PetscReal time, const PetscReal *c,
>> const PetscReal *n, const PetscScalar *a_xI, const PetscScalar *a_xG, void
>> *ctx)
>> {
>> PetscObjectUseFortranCallback((PetscDS)ctx, bocofunc,
>> (PetscReal*,const PetscReal*,const
>> PetscReal*,const PetscScalar*,const PetscScalar*,void*,PetscErrorCode*),
>> (&time,c,n,a_xI,a_xG,ctx,&ierr));
>> }
>> static PetscErrorCode ourbocofunc_time(PetscReal time, const PetscReal *c,
>> const PetscReal *n, const PetscScalar *a_xI, const PetscScalar *a_xG, void
>> *ctx)
>> {
>> PetscObjectUseFortranCallback((PetscDS)ctx, bocofunc_time,
>> (PetscReal*,const PetscReal*,const
>> PetscReal*,const PetscScalar*,const PetscScalar*,void*,PetscErrorCode*),
>> (&time,c,n,a_xI,a_xG,ctx,&ierr));
>> }
>> PETSC_EXTERN void petscdsaddboundary_(PetscDS *prob,
>> DMBoundaryConditionType *type, char *name, char *labelname, PetscInt
>> *field, PetscInt *numcomps, PetscInt *comps,
>> void (*bcFunc)(void),
>> void (*bcFunc_t)(void),
>> PetscInt *numids, const PetscInt
>> *ids, void *ctx, PetscErrorCode *ierr,
>> PETSC_FORTRAN_CHARLEN_T namelen,
>> PETSC_FORTRAN_CHARLEN_T labelnamelen)
>> {
>> char *newname, *newlabelname;
>> FIXCHAR(name, namelen, newname);
>> FIXCHAR(labelname, labelnamelen, newlabelname);
>> *ierr = PetscObjectSetFortranCallback((PetscObject)*prob,
>> PETSC_FORTRAN_CALLBACK_CLASS, &bocofunc, (PetscVoidFunction)bcFunc,
>> ctx);
>> *ierr = PetscObjectSetFortranCallback((PetscObject)*prob,
>> PETSC_FORTRAN_CALLBACK_CLASS, &bocofunc_time, (PetscVoidFunction)bcFunc_t,
>> ctx);
>> *ierr = PetscDSAddBoundary(*prob, *type, newname, newlabelname,
>> *field, *numcomps, comps,
>> (void (*)(void))ourbocofunc,
>> (void (*)(void))ourbocofunc_time,
>> *numids, ids, *prob);
>> FREECHAR(name, newname);
>> FREECHAR(labelname, newlabelname);
>> }
>>
>>
>>
>> but when I call it in the program, with adequate routines, I obtain the
>> following error :
>>
>> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------[0]PETSC ERROR: Corrupt argument: https://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC ERROR: Fortran callback not set on this object[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.[0]PETSC ERROR: Petsc Development GIT revision: v3.14.2-297-gf36a7edeb8 GIT Date: 2020-12-18 04:42:53 +0000[0]PETSC ERROR: ../../../bin/eulerian3D on a named macbook-pro-de-thibault.home by tbridel Sun Dec 20 15:05:15 2020[0]PETSC ERROR: Configure options --with-clean=0 --prefix=/Users/tbridel/Documents/1-CODES/04-PETSC/build --with-make-np=2 --with-windows-graphics=0 --with-debugging=0 --download-fblaslapack --download-mpich-shared=0 --with-x=0 --with-pthread=0 --with-valgrind=0 --PETSC_ARCH=macosx --with-fc=/usr/local/Cellar/open-mpi/4.0.1_2/bin/mpifort --with-cc=/usr/local/Cellar/open-mpi/4.0.1_2/bin/mpicc --with-cxx=/usr/local/Cellar/open-mpi/4.0.1_2/bin/mpic++ --with-openmp=0 --download-hypre=yes --download-sowing=yes --download-metis=yes --download-parmetis=yes --download-triangle=yes --download-tetgen=yes --download-ctetgen=yes --download-p4est=yes --download-zlib=yes --download-c2html=yes --download-eigen=yes --download-pragmatic=yes --with-hdf5-dir=/usr/local/Cellar/hdf5/1.10.5_1 --with-cmake-dir=/usr/local/Cellar/cmake/3.15.3[0]PETSC ERROR: #1 PetscObjectGetFortranCallback() line 258 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/sys/objects/inherit.c[0]PETSC ERROR: #2 ourbocofunc() line 141 in /Users/tbridel/Documents/1-CODES/59-EULERIAN3D/sources/petsc_wrapping/wrapper_petsc.c[0]PETSC ERROR: #3 DMPlexInsertBoundaryValuesRiemann() line 989 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/dm/impls/plex/plexfem.c[0]PETSC ERROR: #4 DMPlexInsertBoundaryValues_Plex() line 1052 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/dm/impls/plex/plexfem.c[0]PETSC ERROR: #5 DMPlexInsertBoundaryValues() line 1142 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/dm/impls/plex/plexfem.c[0]PETSC ERROR: #6 DMPlexComputeResidual_Internal() line 4524 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/dm/impls/plex/plexfem.c[0]PETSC ERROR: #7 DMPlexTSComputeRHSFunctionFVM() line 74 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/ts/utils/dmplexts.c[0]PETSC ERROR: #8 ourdmtsrhsfunc() line 186 in /Users/tbridel/Documents/1-CODES/59-EULERIAN3D/sources/petsc_wrapping/wrapper_petsc.c[0]PETSC ERROR: #9 TSComputeRHSFunction_DMLocal() line 105 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/ts/utils/dmlocalts.c[0]PETSC ERROR: #10 TSComputeRHSFunction() line 653 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/ts/interface/ts.c[0]PETSC ERROR: #11 TSSSPStep_RK_3() line 120 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/ts/impls/explicit/ssp/ssp.c[0]PETSC ERROR: #12 TSStep_SSP() line 208 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/ts/impls/explicit/ssp/ssp.c[0]PETSC ERROR: #13 TSStep() line 3757 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/ts/interface/ts.c[0]PETSC ERROR: #14 TSSolve() line 4154 in /Users/tbridel/Documents/1-CODES/04-PETSC/src/ts/interface/ts.c[0]PETSC ERROR: #15 User provided function() line 0 in User file
>>
>>
>> Second is about the DMProjectFunction wrapper, that I have done so :
>>
>> static PetscErrorCode ourdmprojfunc(PetscInt dim, PetscReal time,
>> PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
>> {
>> PetscObjectUseFortranCallback((DM)ctx, dmprojfunc,
>>
>> (PetscInt*,PetscReal*,PetscReal*,PetscInt*,PetscScalar*,void*,PetscErrorCode*),
>> (&dim,&time,x,&Nf,u,_ctx,&ierr))
>> }
>> PETSC_EXTERN void dmprojectfunction_(DM *dm, PetscReal *time,
>> void
>> (*func)(PetscInt*,PetscReal*,PetscReal*,PetscInt*,PetscScalar*,void*,PetscErrorCode*),
>> void *ctx, InsertMode *mode, Vec X,
>> PetscErrorCode *ierr)
>> {
>> PetscErrorCode (*funcarr[1]) (PetscInt dim, PetscReal time, PetscReal
>> x[], PetscInt Nf, PetscScalar *u, void *ctx);
>> *ierr = PetscObjectSetFortranCallback((PetscObject)*dm,
>> PETSC_FORTRAN_CALLBACK_CLASS, &dmprojfunc, (PetscVoidFunction)func, ctx);
>> funcarr[0] = ourdmprojfunc;
>> *ierr = DMProjectFunction(*dm, *time, funcarr, &ctx, *mode, X);
>> }
>>
>>
>> This time there is no error because I cannot reach this point in the
>> program, but I am not sure anyways how to write this wrapper, especially
>> because of the double pointers that DMProjectFunction takes as arguments.
>>
>> Does anyone have any idea what could be going wrong with those two
>> wrappers ?
>>
>> Thank you very much in advance !!
>>
>> Thibault
>>
>> Le ven. 18 déc. 2020 à 11:02, Thibault Bridel-Bertomeu <
>> thibault.bridelbertomeu at gmail.com> a écrit :
>>
>>> Aah that is a nice trick, I was getting ready to fork, clone the fork and
>>> redo the work, but that worked fine ! Thank you Barry !
>>>
>>> The MR will appear in a little while !
>>>
>>> Thibault
>>>
>>>
>>> Le ven. 18 déc. 2020 à 10:16, Barry Smith <bsmith at petsc.dev> a écrit :
>>>
>>>>
>>>> Good question. There is a trick to limit the amount of work you need
>>>> to do with a new fork after you have already made changes with a PETSc
>>>> clone, but it looks like we do not document this clearly in the webpages.
>>>> (I couldn't find it).
>>>>
>>>> Yes, you do need to make a fork, but after you have made the fork on
>>>> the GitLab website (and have done nothing on your machine) edit the file
>>>> $PETSC_DIR/.git/config for your clone on your machine
>>>>
>>>> Locate the line that has url = git at gitlab.com:petsc/petsc.git (this
>>>> may have an https at the beginning of the line)
>>>>
>>>> Change this line to point to the fork url instead with git@ not
>>>> https, which will be pretty much the same URL but with your user id instead
>>>> of petsc in the address. Then git push and it will push to your fork.
>>>>
>>>> Now you changes will be in your fork and you can make the MR from your
>>>> fork URL on Gitlab. (In other words this editing trick converts your PETSc
>>>> clone on your machine to a PETSc fork).
>>>>
>>>> I hope I have explained this clearly enough it goes smoothly.
>>>>
>>>> Barry
>>>>
>>>>
>>>>
>>>> On Dec 18, 2020, at 3:00 AM, Thibault Bridel-Bertomeu <
>>>> thibault.bridelbertomeu at gmail.com> wrote:
>>>>
>>>> Hello Barry,
>>>>
>>>> I'll start the MR as soon as possible then so that specialists can
>>>> indeed have a look. Do I have to fork PETSc to start a MR or are PETSc repo
>>>> settings such that can I push a branch from the PETSc clone I got ?
>>>>
>>>> Thibault
>>>>
>>>>
>>>> Le mer. 16 déc. 2020 à 07:47, Barry Smith <bsmith at petsc.dev> a écrit :
>>>>
>>>>>
>>>>> Thibault,
>>>>>
>>>>> A subdirectory for the example is fine; we have other examples that
>>>>> use subdirectories and multiple files.
>>>>>
>>>>> Note: even if you don't have something completely working you can
>>>>> still make MR and list it as DRAFT request for comments, some other PETSc
>>>>> members who understand the packages you are using and Fortran better than I
>>>>> may be able to help as you develop the code.
>>>>>
>>>>> Barry
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On Dec 16, 2020, at 12:35 AM, Thibault Bridel-Bertomeu <
>>>>> thibault.bridelbertomeu at gmail.com> wrote:
>>>>>
>>>>> Hello everyone,
>>>>>
>>>>> Thank you Barry for the feedback.
>>>>> OK, yes I'll work up an MR as soon as I have got something working. By
>>>>> the way, does the fortran-version of the example have to be a single file ?
>>>>> If my push contains a directory with several files (different modules and
>>>>> the main), and the Makefile that goes with it, is that ok ?
>>>>>
>>>>> Thibault Bridel-Bertomeu
>>>>>
>>>>>
>>>>> Le mer. 16 déc. 2020 à 04:46, Barry Smith <bsmith at petsc.dev> a écrit :
>>>>>
>>>>>>
>>>>>> This is great. If you make a branch off of the PETSc git repository
>>>>>> with these additions and work on ex11 you can make a merge request and we
>>>>>> can run the code easily on all our test systems (for security reasons one
>>>>>> of use needs to launch the tests from your MR).
>>>>>> https://docs.petsc.org/en/latest/developers/integration/
>>>>>>
>>>>>> Barry
>>>>>>
>>>>>>
>>>>>> On Dec 15, 2020, at 5:35 AM, Thibault Bridel-Bertomeu <
>>>>>> thibault.bridelbertomeu at gmail.com> wrote:
>>>>>>
>>>>>> Hello everyone,
>>>>>>
>>>>>> So far, I have the wrappers in the files attached to this e-mail. I
>>>>>> still do not know if they work properly - at least the code compiles and
>>>>>> the calls to the wrapped-subroutine do not fail - but I wanted to put this
>>>>>> here in case someone sees something really wrong with it already.
>>>>>>
>>>>>> Thank you again for your help, I'll try to post updates of the F90
>>>>>> version of ex11 regularly in this thread.
>>>>>>
>>>>>> Stay safe,
>>>>>>
>>>>>> Thibault Bridel-Bertomeu
>>>>>>
>>>>>> Le dim. 13 déc. 2020 à 16:39, Jed Brown <jed at jedbrown.org> a écrit :
>>>>>>
>>>>>>> Thibault Bridel-Bertomeu <thibault.bridelbertomeu at gmail.com> writes:
>>>>>>>
>>>>>>> > Thank you Mark for your answer.
>>>>>>> >
>>>>>>> > I am not sure what you think could be in the setBC1 routine ? How
>>>>>>> to make
>>>>>>> > the connection with the PetscDS ?
>>>>>>> >
>>>>>>> > On the other hand, I actually found after a while TSMonitorSet has a
>>>>>>> > fortran wrapper, and it does take as arguments two function
>>>>>>> pointers, so I
>>>>>>> > guess it is possible ? Although I am not sure exactly how to play
>>>>>>> with the
>>>>>>> > PetscObjectSetFortranCallback & PetscObjectUseFortranCallback
>>>>>>> macros -
>>>>>>> > could anybody advise please ?
>>>>>>>
>>>>>>> tsmonitorset_ is a good example to follow. In your file, create one
>>>>>>> of these static structs with a member for each callback. These are IDs that
>>>>>>> will be used as keys for Fortran callbacks and their contexts. The salient
>>>>>>> parts of the file are below.
>>>>>>>
>>>>>>> static struct {
>>>>>>> PetscFortranCallbackId prestep;
>>>>>>> PetscFortranCallbackId poststep;
>>>>>>> PetscFortranCallbackId rhsfunction;
>>>>>>> PetscFortranCallbackId rhsjacobian;
>>>>>>> PetscFortranCallbackId ifunction;
>>>>>>> PetscFortranCallbackId ijacobian;
>>>>>>> PetscFortranCallbackId monitor;
>>>>>>> PetscFortranCallbackId mondestroy;
>>>>>>> PetscFortranCallbackId transform;
>>>>>>> #if defined(PETSC_HAVE_F90_2PTR_ARG)
>>>>>>> PetscFortranCallbackId function_pgiptr;
>>>>>>> #endif
>>>>>>> } _cb;
>>>>>>>
>>>>>>> /*
>>>>>>> Note ctx is the same as ts so we need to get the Fortran context
>>>>>>> out of the TS; this gets put in _ctx using the callback ID
>>>>>>> */
>>>>>>> static PetscErrorCode ourmonitor(TS ts,PetscInt i,PetscReal d,Vec
>>>>>>> v,void *ctx)
>>>>>>> {
>>>>>>>
>>>>>>> PetscObjectUseFortranCallback(ts,_cb.monitor,(TS*,PetscInt*,PetscReal*,Vec
>>>>>>> *,void*,PetscErrorCode*),(&ts,&i,&d,&v,_ctx,&ierr));
>>>>>>> }
>>>>>>>
>>>>>>> Then follow as in tsmonitorset_, which sets two callbacks.
>>>>>>>
>>>>>>> PETSC_EXTERN void tsmonitorset_(TS *ts,void
>>>>>>> (*func)(TS*,PetscInt*,PetscReal*,Vec*,void*,PetscErrorCode*),void
>>>>>>> *mctx,void (*d)(void*,PetscErrorCode*),PetscErrorCode *ierr)
>>>>>>> {
>>>>>>> CHKFORTRANNULLFUNCTION(d);
>>>>>>> if ((PetscVoidFunction)func == (PetscVoidFunction)
>>>>>>> tsmonitordefault_) {
>>>>>>> *ierr = TSMonitorSet(*ts,(PetscErrorCode
>>>>>>> (*)(TS,PetscInt,PetscReal,Vec,void*))TSMonitorDefault,*(PetscViewerAndFormat**)mctx,(PetscErrorCode
>>>>>>> (*)(void **))PetscViewerAndFormatDestroy);
>>>>>>> } else {
>>>>>>> *ierr =
>>>>>>> PetscObjectSetFortranCallback((PetscObject)*ts,PETSC_FORTRAN_CALLBACK_CLASS,&_cb.monitor,(PetscVoidFunction)func,mctx);
>>>>>>> *ierr =
>>>>>>> PetscObjectSetFortranCallback((PetscObject)*ts,PETSC_FORTRAN_CALLBACK_CLASS,&_cb.mondestroy,(PetscVoidFunction)d,mctx);
>>>>>>> *ierr = TSMonitorSet(*ts,ourmonitor,*ts,ourmonitordestroy);
>>>>>>> }
>>>>>>> }
>>>>>>>
>>>>>> <wrapper_petsc.h90><wrapper_petsc.c>
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>
>>
More information about the petsc-users
mailing list