[petsc-dev] EM solver in XGC
Mark Adams
mfadams at lbl.gov
Tue Apr 14 18:39:57 CDT 2015
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 <knepley at gmail.com> wrote:
> On Tue, Apr 14, 2015 at 5:19 PM, Mark Adams <mfadams at lbl.gov> 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 <mfadams at lbl.gov> 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 <mfadams at lbl.gov> 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 <mfadams at lbl.gov> 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 <knepley at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> On Tue, Apr 14, 2015 at 12:40 PM, Matthew Knepley <knepley at gmail.com>
>>>>>> wrote:
>>>>>>
>>>>>>> On Tue, Apr 14, 2015 at 12:33 PM, Mark Adams <mfadams at lbl.gov>
>>>>>>> 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 <mfadams at lbl.gov>
>>>>>>>> wrote:
>>>>>>>>
>>>>>>>>> OK,
>>>>>>>>>
>>>>>>>>> On Tue, Apr 14, 2015 at 1:25 PM, Matthew Knepley <
>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>
>>>>>>>>>> On Tue, Apr 14, 2015 at 12:23 PM, Mark Adams <mfadams at lbl.gov>
>>>>>>>>>> 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 <
>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>
>>>>>>>>>>>> On Tue, Apr 14, 2015 at 12:13 PM, Mark Adams <mfadams at lbl.gov>
>>>>>>>>>>>> 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 <
>>>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>>>
>>>>>>>>>>>>>> On Tue, Apr 14, 2015 at 11:16 AM, Mark Adams <mfadams at lbl.gov
>>>>>>>>>>>>>> > 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 <
>>>>>>>>>>>>>>> mfadams at lbl.gov> 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 <
>>>>>>>>>>>>>>>> mfadams at lbl.gov> 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 <
>>>>>>>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>> On Tue, Apr 14, 2015 at 10:35 AM, Mark Adams <
>>>>>>>>>>>>>>>>>> mfadams at lbl.gov> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20150414/bd7f7478/attachment.html>
More information about the petsc-dev
mailing list