Hi Jed, Thanks for your message. I implemented everything in C as you suggested and it works fine except for one thing : the ts->vec_sol does not seem to get updated when seen from the C code (it is on the Fortran side though the solution is correct). As a result, the time step (that uses among other things the max velocity in the domain) is always at the value it gets from the initial solution. Any idea why ts->vec_sol does not seem to be updated ? (I checked the stepnum and time is updated though , when accessed with TSGetTime for instance).
Cheers, Thibault Le lun. 28 déc. 2020 à 15:30, Jed Brown <[email protected]> a écrit : > Thibault Bridel-Bertomeu <[email protected]> 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 <[email protected]> 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 < > >> [email protected]> 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 < > >> [email protected]> 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 <[email protected]> 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 = [email protected]: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 < > >>>> [email protected]> 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 <[email protected]> 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 < > >>>>> [email protected]> 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 <[email protected]> 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 < > >>>>>> [email protected]> 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 <[email protected]> a écrit > : > >>>>>> > >>>>>>> Thibault Bridel-Bertomeu <[email protected]> > 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> > >>>>>> > >>>>>> > >>>>>> > >>>>> > >>>> > >> > -- Thibault Bridel-Bertomeu — Eng, MSc, PhD Research Engineer CEA/CESTA 33114 LE BARP Tel.: (+33)557046924 Mob.: (+33)611025322 Mail: [email protected]
