[petsc-dev] EM solver in XGC
Mark Adams
mfadams at lbl.gov
Tue Apr 14 13:54:49 CDT 2015
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
>>>
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20150414/3c33dcb1/attachment.html>
More information about the petsc-dev
mailing list