OK, I use: call MatConvert(a_ts%FJacobian,MATAIJ,MAT_INITIAL_MATRIX,a_ts%FJacobian2,ierr);CHKERRQ(ierr)
I added code to src/snes/examples/tutorials/ex73f90t.F90 that reproduces this error, in next. I'll look at it. Mark On Tue, Apr 14, 2015 at 6:33 PM, Matthew Knepley <[email protected]> wrote: > On Tue, Apr 14, 2015 at 5:19 PM, Mark Adams <[email protected]> wrote: > >> Matt, in looking at MatConvert_Nest_AIJ, the code looks like it is 2 >> years old. >> >> I need a MatConvert_Nest_MPIAIJ, should I write this or does it exist >> someplace? >> > > This is MPIAIJ. Look at the code. > > Matt > > >> >> Mark >> >> >> >> On Tue, Apr 14, 2015 at 2:54 PM, Mark Adams <[email protected]> wrote: >> >>> Also, this is with MPI so I need MPIAIJ. Is that supposed to work? >>> >>> On Tue, Apr 14, 2015 at 2:53 PM, Mark Adams <[email protected]> wrote: >>> >>>> It looks like MatConvert is looking for "MatConvert_nest_aij_C" >>>> >>>> but I see: >>>> >>>> > git grep MatConvert_Nest >>>> impls/nest/matnest.c:#define __FUNCT__ "MatConvert_Nest_AIJ" >>>> impls/nest/matnest.c:PETSC_EXTERN PetscErrorCode >>>> MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) >>>> >>>> Is there a problems with "Nest" vs "nest" here? >>>> >>>> Mark >>>> >>>> >>>> On Tue, Apr 14, 2015 at 2:46 PM, Mark Adams <[email protected]> wrote: >>>> >>>>> I have tried >>>>> >>>>> call >>>>> MatConvert(a_ts%FJacobian,MATMPIAIJ,MAT_INITIAL_MATRIX,a_ts%FJacobian2,ierr);CHKERRQ(ierr) >>>>> >>>>> and >>>>> >>>>> call >>>>> MatConvert(a_ts%FJacobian,MATAIJ,MAT_INITIAL_MATRIX,a_ts%FJacobian2,ierr);CHKERRQ(ierr) >>>>> >>>>> in 'next' and it is not finding the specialized routine. >>>>> >>>>> I will look into it. >>>>> >>>>> Mark >>>>> >>>>> >>>>> On Tue, Apr 14, 2015 at 1:41 PM, Matthew Knepley <[email protected]> >>>>> wrote: >>>>> >>>>>> On Tue, Apr 14, 2015 at 12:40 PM, Matthew Knepley <[email protected]> >>>>>> wrote: >>>>>> >>>>>>> On Tue, Apr 14, 2015 at 12:33 PM, Mark Adams <[email protected]> >>>>>>> wrote: >>>>>>> >>>>>>>> MatNest does not support convert because it does not have a getrow. >>>>>>>> >>>>>>> >>>>>>> I think this got fixed in dev. Anyhow, this supports my contention >>>>>>> that MatNest is a piece of shit. >>>>>>> >>>>>> >>>>>> Use 'next', it has >>>>>> >>>>>> PETSC_EXTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType >>>>>> newtype,MatReuse reuse,Mat *newmat) >>>>>> >>>>>> so it should work. >>>>>> >>>>>> Matt >>>>>> >>>>>> >>>>>>> >>>>>>> >>>>>>>> Should we use a 4 block fieldsplit? >>>>>>>> >>>>>>> >>>>>>> We give up our ability to check using exact solves, but go ahead and >>>>>>> try it if you think it will make >>>>>>> the physicists happy for now. >>>>>>> >>>>>>> Thanks, >>>>>>> >>>>>>> Matt >>>>>>> >>>>>>> >>>>>>>> On Tue, Apr 14, 2015 at 1:26 PM, Mark Adams <[email protected]> >>>>>>>> wrote: >>>>>>>> >>>>>>>>> OK, >>>>>>>>> >>>>>>>>> On Tue, Apr 14, 2015 at 1:25 PM, Matthew Knepley < >>>>>>>>> [email protected]> wrote: >>>>>>>>> >>>>>>>>>> On Tue, Apr 14, 2015 at 12:23 PM, Mark Adams <[email protected]> >>>>>>>>>> wrote: >>>>>>>>>> >>>>>>>>>>> I moved this to after the setup, which it has to do because it >>>>>>>>>>> does not have the prefix or set options. >>>>>>>>>>> >>>>>>>>>>> It is not failing with not finding the IS. Will MatNest be able >>>>>>>>>>> to merge two matrices together like this? >>>>>>>>>>> >>>>>>>>>>> I create a 4x4 MatNext but have just two fields in fieldSplit. >>>>>>>>>>> >>>>>>>>>> >>>>>>>>>> Yep, I believe that will work. If not, just convert the MatNest >>>>>>>>>> to an AIJ for now with MatConvert(). >>>>>>>>>> >>>>>>>>>> Matt >>>>>>>>>> >>>>>>>>>> >>>>>>>>>>> On Tue, Apr 14, 2015 at 1:18 PM, Matthew Knepley < >>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>> >>>>>>>>>>>> On Tue, Apr 14, 2015 at 12:13 PM, Mark Adams <[email protected]> >>>>>>>>>>>> wrote: >>>>>>>>>>>> >>>>>>>>>>>>> I get an error where it looks like the number fields in field >>>>>>>>>>>>> split is < 2. Do I need to set the size? >>>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> No, each call to SetIS() just appends to the field list. >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>>> I use this code but I do this before I have setup the TS. >>>>>>>>>>>>> Should I do this after I setup the TS? >>>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> This is why I hate using the API since it is difficult to get >>>>>>>>>>>> all the setup done at the right time. First, >>>>>>>>>>>> when you call PCFieldSplitSetIS(), is the PC type set? >>>>>>>>>>>> >>>>>>>>>>>> Thanks, >>>>>>>>>>>> >>>>>>>>>>>> Matt >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>>> ! setup field split >>>>>>>>>>>>> call TSGetSNES(a_ts,snes,ierr);CHKERRQ(ierr) >>>>>>>>>>>>> call SNESGetKSP(snes,ksp,ierr);CHKERRQ(ierr) >>>>>>>>>>>>> call KSPGetPC(ksp,pc,ierr);CHKERRQ(ierr) >>>>>>>>>>>>> >>>>>>>>>>>>> call >>>>>>>>>>>>> MatGetSize(a_ts%mass,PETSC_NULL_INTEGER,npetsc,ierr);CHKERRQ(ierr) >>>>>>>>>>>>> npetsc = 2*npetsc ! two blocks in each partition >>>>>>>>>>>>> j = 0 >>>>>>>>>>>>> call >>>>>>>>>>>>> ISCreateStride(w_comm,npetsc,j,ione,is,ierr);CHKERRQ(ierr) >>>>>>>>>>>>> call PCFieldSplitSetIS(pc,'dt',is,ierr);CHKERRQ(ierr) >>>>>>>>>>>>> call ISDestroy(is,ierr);CHKERRQ(ierr) >>>>>>>>>>>>> j = npetsc >>>>>>>>>>>>> call >>>>>>>>>>>>> ISCreateStride(w_comm,npetsc,j,ione,is,ierr);CHKERRQ(ierr) >>>>>>>>>>>>> call PCFieldSplitSetIS(pc,'aux',is,ierr);CHKERRQ(ierr) >>>>>>>>>>>>> call ISDestroy(is,ierr);CHKERRQ(ierr) >>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> On Tue, Apr 14, 2015 at 12:24 PM, Matthew Knepley < >>>>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>>>> >>>>>>>>>>>>>> On Tue, Apr 14, 2015 at 11:16 AM, Mark Adams <[email protected] >>>>>>>>>>>>>> > wrote: >>>>>>>>>>>>>> >>>>>>>>>>>>>>> Should I just use >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y) >>>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> Yep, that is exactly what I use it for. >>>>>>>>>>>>>> >>>>>>>>>>>>>> Matt >>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>>> ? >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> On Tue, Apr 14, 2015 at 12:15 PM, Mark Adams < >>>>>>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> I am thinking that I should create an IS that is my local >>>>>>>>>>>>>>>> part of each field and then get a "sub vector" with that >>>>>>>>>>>>>>>> somehow. Looking >>>>>>>>>>>>>>>> into this now. >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> On Tue, Apr 14, 2015 at 12:00 PM, Mark Adams < >>>>>>>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> Whoop, old code still there. >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> I now have to get access to vectors somehow. >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> I now get vectors for a field and do stuff with them like: >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> subroutine >>>>>>>>>>>>>>>>> FormIFunction(ts,t_dummy,a_XX,a_Xdot,a_FF,a_ts,ierr) >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> call TSGetDM(ts,dm,ierr);CHKERRQ(ierr) >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> call >>>>>>>>>>>>>>>>> DMCompositeGetAccessArray(dm,a_FF,ifour,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr) >>>>>>>>>>>>>>>>> call >>>>>>>>>>>>>>>>> DMCompositeGetAccessArray(dm,a_Xdot,ifour,PETSC_NULL_INTEGER,Xdotsub,ierr);CHKERRQ(ierr) >>>>>>>>>>>>>>>>> one = 1d0 >>>>>>>>>>>>>>>>> call VecAXPY(Fsub(1),one,Xdotsub(1),ierr);CHKERRQ(ierr) >>>>>>>>>>>>>>>>> call VecAXPY(Fsub(2),one,Xdotsub(2),ierr);CHKERRQ(ierr) >>>>>>>>>>>>>>>>> call >>>>>>>>>>>>>>>>> DMCompositeRestoreAccessArray(dm,a_Xdot,ifour,PETSC_NULL_INTEGER,Xdotsub,ierr);CHKERRQ(ierr) >>>>>>>>>>>>>>>>> call >>>>>>>>>>>>>>>>> DMCompositeRestoreAccessArray(dm,a_FF,ifour,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr) >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> call >>>>>>>>>>>>>>>>> MatMultAdd(a_ts%FJacobian,a_XX,a_FF,a_FF,ierr);CHKERRQ(ierr) >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> I need to get rid of DMCompositeGetAccessArray. What do I >>>>>>>>>>>>>>>>> do? >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> Thanks, >>>>>>>>>>>>>>>>> Mark >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> On Tue, Apr 14, 2015 at 11:36 AM, Matthew Knepley < >>>>>>>>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> On Tue, Apr 14, 2015 at 10:35 AM, Mark Adams < >>>>>>>>>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> And is this the solver (appended). >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> Do not call SNESSetDM(). We are not using DM at all in >>>>>>>>>>>>>>>>>> this code. >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> Matt >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> I am getting an error at this last line: >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() >>>>>>>>>>>>>>>>>>> @*/ >>>>>>>>>>>>>>>>>>> PetscErrorCode SNESSetDM(SNES snes,DM dm) >>>>>>>>>>>>>>>>>>> { >>>>>>>>>>>>>>>>>>> PetscErrorCode ierr; >>>>>>>>>>>>>>>>>>> KSP ksp; >>>>>>>>>>>>>>>>>>> DMSNES sdm; >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> PetscFunctionBegin; >>>>>>>>>>>>>>>>>>> PetscValidHeaderSpecific(snes,SNES_CLASSID,1); >>>>>>>>>>>>>>>>>>> if (dm) {ierr = >>>>>>>>>>>>>>>>>>> PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} >>>>>>>>>>>>>>>>>>> if (snes->dm) { /* Move the DMSNES >>>>>>>>>>>>>>>>>>> context over to the new DM unless the new DM already has >>>>>>>>>>>>>>>>>>> one */ >>>>>>>>>>>>>>>>>>> if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) >>>>>>>>>>>>>>>>>>> { >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> The problem is that the argument 'dm' to this method is >>>>>>>>>>>>>>>>>>> NULL so I get a segv at the last term of this conditional. >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> Mark >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> -adv_ts_type beuler >>>>>>>>>>>>>>>>>>> -adv_ts_monitor >>>>>>>>>>>>>>>>>>> -adv_ts_view >>>>>>>>>>>>>>>>>>> -adv_ts_max_step 1 >>>>>>>>>>>>>>>>>>> -adv_snes_max_it 1 # make >> 1 for >>>>>>>>>>>>>>>>>>> nonlinear solve >>>>>>>>>>>>>>>>>>> -adv_snes_lag_preconditioner -1 # make > 0 to >>>>>>>>>>>>>>>>>>> update PC >>>>>>>>>>>>>>>>>>> -adv_snes_lag_jacobian -1 # make > 0 for >>>>>>>>>>>>>>>>>>> nonlinear Jacobian >>>>>>>>>>>>>>>>>>> -adv_snes_monitor >>>>>>>>>>>>>>>>>>> -adv_snes_converged_reason >>>>>>>>>>>>>>>>>>> -adv_ksp_monitor >>>>>>>>>>>>>>>>>>> -adv_ksp_converged_reason >>>>>>>>>>>>>>>>>>> -adv_ksp_type richardson >>>>>>>>>>>>>>>>>>> -adv_pc_type fieldsplit >>>>>>>>>>>>>>>>>>> #-adv_pc_fieldsplit_type multiplicative >>>>>>>>>>>>>>>>>>> -adv_pc_fieldsplit_type schur >>>>>>>>>>>>>>>>>>> -adv_pc_fieldsplit_schur_fact_type full >>>>>>>>>>>>>>>>>>> -adv_pc_fieldsplit_schur_precondition full >>>>>>>>>>>>>>>>>>> #-adv_pc_fieldsplit_0_fields 0,1 >>>>>>>>>>>>>>>>>>> #-adv_pc_fieldsplit_1_fields 2,3 >>>>>>>>>>>>>>>>>>> -adv_fieldsplit_dt_pc_type lu >>>>>>>>>>>>>>>>>>> -adv_fieldsplit_aux_pc_type lu >>>>>>>>>>>>>>>>>>> -adv_fieldsplit_dt_pc_factor_mat_solver_package superlu >>>>>>>>>>>>>>>>>>> -adv_fieldsplit_aux_pc_factor_mat_solver_package superlu >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>> -- >>>>>>>>>>>>>>>>>> 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 >>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> -- >>>>>>>>>>>>>> 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 >>>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> -- >>>>>>>>>>>> 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 >>>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> -- >>>>>>>>>> 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 >>>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>> >>>>>>> >>>>>>> >>>>>>> -- >>>>>>> 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 >>>>>>> >>>>>> >>>>>> >>>>>> >>>>>> -- >>>>>> 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 >>>>>> >>>>> >>>>> >>>> >>> >> > > > -- > 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 >
