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