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 >
