[petsc-dev] EM solver in XGC
Matthew Knepley
knepley at gmail.com
Tue Apr 14 17:33:59 CDT 2015
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/b979cbe6/attachment.html>
More information about the petsc-dev
mailing list